wspr4fskdemod.m

%%
close all
clear all
clc
%load('wspr4fsk.mat'); %Import received FSK data set bpfsk1.mat


%%
wsprcvd=audiorecorder(48000,8,1);

%%

    %zz=120
            ii=1;
            while ii<120
              disp('.')
              pause(1)
            
                clk=clock;

                %if clk(6)<=1&&mod(clk(5),2)==0
                if clk(6)>=57&&mod(clk(5),2)==1
                disp('Rx!')
                recordblocking(wsprcvd,115)
                ii=120;
                elseif mod(clk(5),2)==1
                    ii=ii+1
                clk = clock;
                disp('waiting to receive')
                end
            end
            
            clip=getaudiodata(wsprcvd);
%% Trim audio clip to correct length

            [xr,locr]=findpeaks(clip);
            meanp = mean(xr);
            [row,colmn]=find(xr>(0.5));
            begin=locr(row(end))-5308416
            clip=clip(begin:locr(row(end)));
            
            plot(clip)
%%
pad=(ones(1,5308416-length(clip)))*mean(clip);
clip=clip';
clip=[pad clip];

f1=1500*(32768/48000)
f2=(1500+(12000/8192))*(32768/48000)
f3=(1500+(24000/8192))*(32768/48000)
f4=(1500+(36000/8192))*(32768/48000)
FSK=clip;
L=162
FS=length(FSK)/L
%FS=1001
%param=[0.0131,0.0133];
param=[0.5,0.51];
Wn1=f1*param/(FS/2);
WN=[0.03 0.04 ]
%BP1=fir1(121,Wn1,'bandpass'); 
BP1=fir1(121,WN,'bandpass'); 
%Wn2=f2*param/(FS/2);
%BP2=fir1(121,Wn2,'bandpass');  
%Wn3=f3*param/(FS/2);
%BP3=fir1(121,Wn3,'bandpass'); 
%Wn4=f4*param/(FS/2);
%BP4=fir1(121,Wn4,'bandpass');  
%secondly,let  received  signal  FSK  passes  the  two  band  pass  filters respectively.
x1=filter(BP1,1,FSK);
%x2=filter(BP2,1,FSK);
%x3=filter(BP3,1,FSK);
%x4=filter(BP4,1,FSK);
%2. Envelope detection %%%%%%%%%%%32
%Mathematically the envelope e(t) of a signal x(t) is defined as the magnitude 
%of the analytic signal(complex signal)/8/. Firstly modify the signal fromreal to 
%complex with the function of Hilbert ( ).

y1=hilbert(x1);
%y2=hilbert(x2);
%y3=hilbert(x3);
%y4=hilbert(x4);

% based on the mathematic theory described before the absolute values of the 
%complex signal can be obtained with the function abs ( ).Envelop of received 
%signal is detected in this way.

envy1=abs(y1);
%envy2=abs(y2);
%envy3=abs(y3);
%envy4=abs(y4);

%%Smoothing function
%%%%%%%%%%%%

%%remove outliers for each pulse
envyp=envy1(1:32768);
mnp=mean(envyp)
for kk=1:32768

	if envy1(kk)>2.4*mnp
	envy1(kk)=mnp;
	end
end

for jj=2:L
    start=((jj-1)*32768+1);
	  fin= (jj*32768);
	envyp=envy1(start:fin);
	mnp=mean(envyp);

	for kk=start:fin

		if envy1(kk)>2.4*mnp
		envy1(kk)=mnp;
		end
	end
end

low=fir1(201,.006,'low');
envy1=filter(low,1,envy1);

%%%%%%%%%%%%%%%%integrate and dump
envyps(1)=sum(envy1(1:32768));


for jj=2:L
    start=((jj-1)*32768+1);
	  fin= (jj*32768);
	envyps(jj)=sum(envy1(start:fin));
end	

maxeps=max(envyps(70:90))/32768
mineps=min(envyps(70:90))/32768

head=maxeps-mineps
headstep=head/3
opts=[mineps (headstep+mineps) (maxeps-headstep) maxeps]


optsround=[]
envyround=size(envy1)
for oo=1:4
optsround(oo)=opts(oo)-envyps(1)/32768
end

est=min(abs(optsround))
for oo=1:4
	if est==abs(optsround(oo))
	%envyround(1:32768)=opts(oo);
    rec4sym(1)=oo-1;
	break
	end	
end



for jj=2:L
    

	for oo=1:4
		optsround(oo)=opts(oo)-envyps(jj-1)/32768;
	end
	
	est=min(abs(optsround))	
	for pp=1:4
	if est==abs(optsround(pp))
	%envyround(jj)=pp-1;
	rec4sym(jj)=4-pp;
	break
	
	end
		pp
			
	end
	
   
end

rec4sym=circshift(rec4sym,-1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Weak Signal Propagation Reporter (WSPR) Recieve Code         %
%      Refernece credit to Dr. Jonathon Y. Cheah (NZ0C)          %
%                   for the use of this code                     %
%           Ultra low baud rate communication study              %
%                                                                %
%               Consult www.wsprnet.org for details              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% clear all Matlab memory buffers
      %close all
      %clear all
      %clc

% get the channel symbol file generated by transmit code
      %signal=textread('signal.dat','%u'); 
      signal=rec4sym';
      signal=vertcat(signal,0,0,0,0,0,0);
      sym=zeros(length(signal),1);
     
% set the initial variable configuration
      MAXBITS=103;        %
      gamma=zeros(MAXBITS-1);           %
      metrics=zeros(4,MAXBITS-1);       %
      tm=zeros(2,MAXBITS-1);         
      nstate=zeros(MAXBITS-1);     %
      
% configure FEC sequence
      npoly1=2^32-221228207;
      npoly2=2^32-463389625;
      
% configure the Fano decoder 
      partab=textread('partab.dat','%u');
      mettab=textread('mettab.dat','%d');
      mettab=reshape(mettab,256,2);
   
% configure the sync sequence
      sync=textread('sync.dat','%u'); 

%   
% Strip the sync sequence from channel signal.
         sym=uint8(sign(signal-sync));
         disp ('wspr signal sync stripped:')
         %sy1=reshape(sym,8,length(sym)/8)';                   %
        % disp (sy1)
         
%  
% De-interleave channel symbols.
%      symbol=deinterleave(sym);
      P=1;
      I=0;
      symbol=zeros(1,168);
      while P < 162
        J=bin2dec(fliplr(bin(fi(I,0,8,0))));
        if J<162
            symbol(P)=sym(J+1);
            P=P+1;
        end
      I=I+1;
      end
      disp ('De-interleaved:')
      sy2=reshape(symbol,8,length(symbol)/8)';                   %
      disp(sy2)
      
%
% Configure fano decoder.
      dat=zeros(11);
      nsym=162;
      nbits=50+31;
      ndelta=50;
      maxcycles=20000;
      ntail=nbits-31;
      i4a=0;
      i4b=0;
      for np =1:nbits
          j=2*np;
          i4a=-symbol(j-1);
          i4b=-symbol(j);
          if (i4a<0)
              i4a=i4a+256;
          end
          if (i4b<0)
              i4b=i4b+256;
          end
          metrics(1,np) = mettab(i4a+1,1) + mettab(i4b+1,1);
          metrics(2,np) = mettab(i4a+1,1) + mettab(i4b+1,2);
          metrics(3,np) = mettab(i4a+1,2) + mettab(i4b+1,1);
          metrics(4,np) = mettab(i4a+1,2) + mettab(i4b+1,2);
      end
      np=1;
      nstate(np)=0;
      
 %     
 % Compute and sort fano branch metrics 
      n=bitand(nstate(np),npoly1);
      n=bitxor(n,bitshift(n,-16));
      m=bitand(bitxor(n,bitshift(n,-8)),255);
      lsym=partab(m+1);
      n=bitand(nstate(np),npoly2);
      n=bitxor(n,bitshift(n,-16));
      m=bitand(bitxor(n,bitshift(n,-8)),255);
      lsym=lsym+lsym+partab(m+1);
      m0=metrics(lsym+1,np);
      m1=metrics(bitxor(3,lsym)+1,np);
  %
  % tm(1,np)=m0 if 0-branch is better
  %         =m1 if 1-branch is better
      if m0>m1
          tm(1,np)= m0;                   
          tm(2,np)=m1;
      else
          tm(1,np)=m1;                    
          tm(2,np)=m0;
          nstate(np)=mod(nstate(np) + 1,2^32) ;     
      end
  
      ii(np)=0;                                
      gamma(np)=0;
      nt=0;
      
  %
  % Start fano with the best branch
      for i=1:nbits*maxcycles-1                
          ngamma=gamma(np) + tm(ii(np)+1,np);
          if ngamma>nt         
              if gamma(np)<(nt+ndelta)
                  nt=nt + ndelta * floor((ngamma-nt)/ndelta);
              end
        
              gamma(np+1)=ngamma;              
              nstate(np+1)=mod(bitshift(nstate(np),1),2^32);
              np=np+1;
              
  %
  % fano decoding done.
              if np == nbits
              break                            
              end
              
  %
  % fano processing.          
              n=mod(bitand(nstate(np),npoly1),2^32);
              n=mod(bitxor(n,bitshift(n,-16)),2^32);
              lsym=partab(bitand(bitxor(n,bitshift(n,-8)),255)+1);
              n=mod(bitand(nstate(np),npoly2),2^32);
              n=mod(bitxor(n,bitshift(n,-16)),2^32);
              lsym=lsym+lsym+partab(bitand(bitxor(n,bitshift(n,-8)),255)+1);
   
   %
   % fano at the tail
              if np >ntail+1
                  tm(1,np)=metrics(lsym+1,np);   
              else
                  m0=metrics(lsym+1,np);
                  m1=metrics(bitxor(3,lsym)+1,np);
                  if m0>m1
                      tm(1,np)=m0;              
                      tm(2,np)=m1;
                  else
                      tm(1,np)=m1;               
                      tm(2,np)=m0;
                      nstate(np)=mod(nstate(np) + 1,2^32);  
                  end
              end
              
  %
  % Start at the best branch
              ii(np)=0  ;                        
               continue     
          end
       npp=0;
       
     while(true)    
       noback=0;
            if np==1
               noback=1;
            end   
            if np>1        
                 if gamma(np-1)<nt
                    noback=1;
                 end
            end
            if(noback)
               nt=nt-ndelta;
                    if ii(np)~=0
                       ii(np)=0;
                       nstate(np)=mod(bitxor(nstate(np),1),32);
                    end
            break   
            end   
            
  % Back up search          
            np=np-1;                            
            if(np<ntail+1 && ii(np)~=1)
  
  %Search next best branch
               ii(np)=ii(np)+1;                   
               nstate(np)=mod(bitxor(nstate(np),1),2^32);
            break
            end
     end
      end
       metric=gamma(np);  
       
  %
  % message code generation 
       nbytes=(nbits+7)/8;                               
       np=8;
for j=1:nbytes-1
    i4a=nstate(np);
      if mod(floor(i4a/128),2)==1
         dat(j) = mod(i4a,128)-128;
        else
         dat(j) = mod(i4a,128);
      end     
    np=np+8;
end

  %
  % Clean-up and display recived message.
   dat(nbytes)=0;
  
   disp ('message code:')
   da=dat(:,1);
   disp(da')
   mess=messageDecode(dat);
   fprintf('Received message = %s\n', mess{:});