%pllcrt: carrier recovery test

clear all

%create pulse shaped received signal

Ts=1/10000;time=2; t=Ts:Ts:time;
M=50; N=length(t)/M;
m=pam(N,4,5);                             % 4-level signal of length N
mup=zeros(1,N*M); mup(1:M:end)=m;         % oversample by M
ps=hamming(M);                            % blip pulse of width M
s=filter(ps,1,mup);                       % convolve pulse shape with data
f0=1000; phoff=-1;                       % carrier freq. and phase
c=cos(2*pi*f0*t+phoff);                   % construct carrier
r=s.*c;                                   % received signal

figure(1)
plotspec(s,Ts)
figure(2)
plotspec(r,Ts)

%send received signal through square and BPF

q=r.^2;                                % square nonlinearity
fl=500;                                % BPF length
yo=4*f0*Ts;
ff=[0 yo-0.01 yo-0.005 yo+0.005 yo+0.01 1];    % BPF center frequency at 2f0
fa=[0 0 1 1 0 0];                     
h=remez(fl,ff,fa);                     % BPF design via remez
r2bpf=filter(h,1,q);                   % filter to give bpf(r^2) 

figure(3)
plotspec(r2bpf,Ts)

%calculating cost function over data set
thta=-3:.1:3;
lent=length(t);
cfl=500;
cff=[0 f0*Ts f0*Ts+0.01 1];
cfa=[1 1 0 0];
ch=remez(cfl,cff,cfa);

for ind=1:length(thta)               % for each possible theta
  x=cos(4*pi*f0*t+2*thta(ind));            % find x with this theta
  gx=filter(ch,1,r2bpf.*x);
  jpll(ind)=sum(gx)/lent;

end

figure(4)
plot(thta,jpll)                  % plot J(theta) vs theta
xlabel('Phase Estimates \theta')
ylabel('Cost Jpll(\theta)' )

%simulating double pll  
fc=f0;mu1=0.001; mu2=0.0002; 
%fc=1.001*f0;mu1=0.005; mu2=0.0005; 
theta=zeros(1,length(t)); theta(1)=0;     % initialize vector for estimates
th2=zeros(1,length(t)); th2(1)=0;  

for k=1:length(t)-1                       % algorithm update
  theta(k+1)=theta(k)-mu1*r2bpf(k)*sin(4*pi*fc*t(k)+2*theta(k));
  th2(k+1)=th2(k)-mu2*r2bpf(k)*sin(4*pi*fc*t(k)+2*theta(k)+2*th2(k));
  end
figure(5)
plot(t,theta)
title('Phase Tracking via the Phase Locked Loop')
xlabel('time'); ylabel('phase offset')
figure(6)
plot(t,th2)
title('Phase Tracking via the Double Phase Locked Loop')
xlabel('time'); ylabel('phase offset')