%KA   4/Feb/2010

%QPSK/16QAM phase tracking
%model for the whole system
%----------------------------------------------------------

clear all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%parameters 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mode = 1;       %1 for qpsk, 0 for 16qam
N  = 5000;      %number of symbols

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%generate random symbols
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if mode == 1 
    alphabet = [+1 -1];
else
    alphabet = [+1 -1 +1/3 -1/3];
end
 
I = randsrc(1,N,alphabet,100);
Q = randsrc(1,N,alphabet,20);

%override for offset marking
I(1000:1025) = 1;  
I(1026:1050) = -1;  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%shaping filter, upsampling x2  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h1 = firrcos(40,.25,.35,1,'rolloff','sqrt');
I_filtered = upfirdn(I,h1,2,1);
Q_filtered = upfirdn(Q,h1,2,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%further upsampling x5 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Iup = resample(I_filtered,5,1);
Qup = resample(Q_filtered,5,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%upconvert to carrier at .2Fs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
carrier = exp(j*2*pi*(0:length(Iup)-1)*0.20000);
IQ = complex(Iup,Qup);
tx_m = real(IQ .* carrier);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%downconvert to baseband 
%leave some carrier residue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
carrier = exp(j*2*pi*(0:length(tx_m)-1)*0.20001);
I_rx1 = tx_m .* real(carrier);
Q_rx1 = tx_m .* -imag(carrier);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%filter off high sideband
%use filtfilt only(ideal)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h2 = fir1(70,.3);
I_rx2 = filtfilt(h2,1,I_rx1); 
Q_rx2 = filtfilt(h2,1,Q_rx1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%downsampling by 5
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I_dn1 = resample(I_rx2,1,5);
Q_dn1 = resample(Q_rx2,1,5);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%downsample by 2, matched filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I_r = upfirdn(I_dn1,h1,1,2);
Q_r = upfirdn(Q_dn1,h1,1,2);

I_r2 = I_r - mean(I_r);        %remove dc
Q_r2 = Q_r - mean(Q_r);

I_r2 = I_r2/max(abs(I_r2));    %scale to unity 
Q_r2 = Q_r2/max(abs(Q_r2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Costas loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

gamma = 1;   %fractional loop filter gain control
alpha = .2;  %integer loop filter gain control

%loop filter coefficients
b = [gamma+alpha -gamma];
a = [1 -1];

for i = 1:N
    
     %slicer
     %%%%%%%%%%%%%%%%%%%%%%%%%
     if mode == 1
  
      if I_r2(i) > 0 
          I_s(i) = 1; 
      else
          I_s(i) = -1; 
      end
      
      if Q_r2(i) > 0 
          Q_s(i) = 1; 
      else
          Q_s(i) = -1; 
      end
   
    else
         %16qam
         if I_r2(i) > 2/3 
             I_s(i) = 1;
         elseif I_r2(i) > 0 
             I_s(i) = 1/3;
         elseif I_r2(i) < -2/3 
             I_s(i) = -1;
         else
             I_s(i) = -1/3;
         end
      
         if Q_r2(i) > 2/3 
             Q_s(i) = 1;
         elseif Q_r2(i) > 0 
             Q_s(i) = 1/3;
         elseif Q_r2(i) < -2/3 
             Q_s(i) = -1;
         else
             Q_s(i) = -1/3;
         end
      
     end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %error detector
    phase_error(i) = Q_r2(i).*I_s(i) - I_r2(i).*Q_s(i);

    %loop filter
    filtered_error(i) = filter(b,a,phase_error(i));

    %accumulator
    accum_error(i) = accum_error(i-1) + filtered_error(i);
    
    %DDS
    carrier = exp(j*2*pi*(0:N-1)* accum_error);
    
    %rotate phase
    
end



      
 




