%ccrt: Costas 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)


%calculating cost function over data set
theta=-3:.1:3;
lent=length(t);
fl=100; ff=[0 .1 .2 1]; fa=[1 1 0 0];
h=remez(fl,ff,fa);                        % LPF design

for ind=1:length(theta)               % for each possible theta
  cth=cos(2*pi*f0*t+theta(ind));
  lpfrcth=filter(h,1,2*r.*cth);
  jc(ind)=0.5*sum(lpfrcth.*lpfrcth)/lent;
end

lpfs=filter(h,1,s);             % cost function scale factor
avsqlpfs=(sum(lpfs.*lpfs))/length(lpfs);  % average squared lpf s
avsqs=(sum(s.*s))/length(s);    % average squared s baseband message signal
figure(3)
plot(theta,jc*2/avsqlpfs)                  % plot J(theta) vs theta
xlabel('Phase Estimates \theta')
ylabel('Cost Jnc(\theta)' )


%simulating Costas
%fc=1.001*f0;mu=0.002; 
fc=f0;mu=0.001;
thc=zeros(1,length(t)); thc(1)=0;     % initialize vector for estimates
zs=zeros(1,fl+1); zc=zeros(1,fl+1);       % initialize buffers for LPFs
for k=1:length(t)-1                       % z's contain past fl+1 inputs
  zs=[2*r(k)*sin(2*pi*fc*t(k)+thc(k)), zs(1:fl)];
  zc=[2*r(k)*cos(2*pi*fc*t(k)+thc(k)), zc(1:fl)];
  lpfs=fliplr(h)*zs'; lpfc=fliplr(h)*zc'; % new output of filters
  thc(k+1)=thc(k)-mu*lpfs*lpfc;       % algorithm update
end
figure(4)
plot(t,thc)
title('Phase Tracking via the Phase Locked Loop')
xlabel('time'); ylabel('phase offset')
