% ddcrt.m: decision directed carrier recovery test

clear all

% pulrecsig.m
Ts=1/10000; time=2; t=Ts:Ts:time;         % time vector
M=50; N=length(t)/M;                      % oversampling factor/amount of data
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

%calculating cost function over data set
thta=-3:.02:3;
lent=length(t);
fl=100; fbe=[0 .2 .3 1]; damps=[1 1 0 0 ]; % parameters for LPF
h=remez(fl,fbe,damps);                     % LPF impulse response
for ind=1:length(thta)               % for each possible theta
  ddx=2*r.*cos(2*pi*f0*t+thta(ind));
  lpddx=filter(h,1,ddx);
  dslpddx=lpddx(0.5*fl+M/2:M:end);
  qlpddx=quantalph(dslpddx,[-3,-1,1,3]);
  jdd(ind)=(qlpddx'-dslpddx)*(qlpddx'-dslpddx)'/length(dslpddx);
end

figure(1)
plot(thta,jdd)                  % plot J(theta) vs theta
xlabel('Phase Estimates \theta')
ylabel('Cost Jdd(\theta)' )

% at the receiver...
% plldd.m: decision directed phase tracking
fzc=zeros(1,fl+1); fzs=zeros(1,fl+1);      % initial state of filters=0
theta=zeros(1,N); theta(1)=-1.27;          % initial phase estimate
mu=.03; j=1; fc=f0;                          % algorithm stepsize mu
for k=1:length(r)
  cc=2*cos(2*pi*fc*t(k)+theta(j));         % cosine for demod
  ss=2*sin(2*pi*fc*t(k)+theta(j));         % sine for demod
  rc=r(k)*cc; rs=r(k)*ss;                  % do the demods                 
  fzc=[rc,fzc(1:fl)]; fzs=[rs,fzs(1:fl)];  % states for LPFs
  x(k)=fliplr(h)*fzc';xder=fliplr(h)*fzs'; % LPFs give x and its derivative
  if mod(0.5*fl+M/2-k,M)==0                % pick correct timing
    qx=quantalph(x(k),[-3,-1,1,3]);        % quantize to nearest symbol
    theta(j+1)=theta(j)-mu*(qx-x(k))*xder; % algorithm update
    j=j+1;
  end
end	

figure(2)
plot(t(mod(0.5*fl+M/2,M):M:end),theta(2:end))
xlabel('time')
ylabel('phase estimates')