clc;
clear all;
close all;

fs=20;     %fixed sample frequency of receiver 
Ts=1/fs;   %fixed sample period of receiver 
fb=3;      %symbol frequency
Tb=1/fb;   %symbol period 
Ti=Tb/2;   %

%% parameter for loop filter
zeta = 0.707; % Damping factor

Wl =0.05*fb; % Lock range defined according to expected input characteristics
Wn = Wl/2/zeta; % Natural frequency calculated from lock range

Kd = 2.27 % Phase detector gain
Ko = 1; % Loop filter gain

C1=((8*zeta*Wn*Tb)/(4+4*zeta*Wn*Tb+Wn*Wn*Tb*Tb))/(Ko*Kd);
C2=((4*Wn*Wn*Ts*Tb)/(4+4*zeta*Wn*Tb+Wn*Wn*Tb*Tb))/(Ko*Kd);


%%

x=importdata('basebandSample.txt');
figure;
plot(x);
title('waveform of baseband sample ');
hold on;
plot(x,'r*');

%%

uk=0; %fractional interval
mkDelta =floor(Ti/Ts+uk ); %increment of basepoint index 

curMk=3;     %current basepoint index 

IdotCnt=1;  %output index for interpolator 
VdCnt=1;    %output index for TED

demodOutCnt=1; %demod decision output index


Discriminator_Out(1)=0; %output of Gardner TED 
PLL_Phase_Part(1)=0;   %gain branch of loop filter
PLL_Freq_Part(1)=0;    %integral branch of loop filter

Control(1)=0;    %error value from loop filter

while curMk <= length(x)-1
   
   y(IdotCnt)=(1-uk)*x(curMk)+uk*x(curMk+1);  %linear interpolate
   
   
   if IdotCnt>= 3 && (mod(IdotCnt ,2)==1)
       
        Discriminator_Out(VdCnt)=y(IdotCnt-1)*[sign(y(IdotCnt))-sign(y(IdotCnt-2))];  %Gardner TED for BPSK
        %Discriminator_Out(VdCnt)=y(IdotCnt-1)*[y(IdotCnt)-y(IdotCnt-2)];
       
        demodOut(demodOutCnt)=y(IdotCnt);
        demodOutCnt=demodOutCnt+1;
        
        %loop filter
       if VdCnt >= 2 

            PLL_Phase_Part(VdCnt)=Discriminator_Out(VdCnt)*C1;
            Control(VdCnt)=PLL_Phase_Part(VdCnt)+PLL_Freq_Part(VdCnt-1);
            PLL_Freq_Part(VdCnt)=Discriminator_Out(VdCnt)*C2+PLL_Freq_Part(VdCnt-1);
           
                      
            %get next increment of basepoint and fractional interval
            mkDelta=floor( Ti/Ts+ Control(VdCnt)+uk );
            uk=mod( Ti/Ts+ Control(VdCnt)+uk ,1); 

       end
       
       VdCnt=VdCnt+1;
       
   end
  
   curMk=curMk+mkDelta; %next basepoint
   IdotCnt=IdotCnt+1; 
 
end


figure;
plot(Discriminator_Out);

figure;
plot(Control);

scatterplot(x);
scatterplot(demodOut(80:length(demodOut)));