% PLL illustration
clear all;
close all;
% define initial phase offset and incoming Cw frequency and Sampling frequency
theta = pi/3;
f=1000;
fs=100000;
% Create the real and imaginary parts of a Cw non-modulated Carrier to be tracked.
k=1:1:1000;
delf=f/20;
Signal=exp(j*(2*pi*k*(f+delf)/fs+theta))+0.01*(rand(1,1000)+j*rand(1,1000));
% initilize PLL Loop
phi_hat(1)=30;
e(1)=0;
phd_output(1)=0;
nco(1)=0
% define Loop filter parameters
kp=0.15; % Proportional constant
ki=0.1; % Integrator constant
% PLL implementation
for n=2:length(Signal)
nco(n)=conj(exp(j*(2*pi*n*f/fs+phi_hat(n-1)))); % Compute nCO
phd_output(n)=imag(Signal(n)*nco(n)); % Complex multiply nCO x input
e(n)=e(n-1)+(kp+ki)*phd_output(n)-ki*phd_output(n-1); % Filter integrator
some(n)=(kp+ki)*phd_output(n)-ki*phd_output(n-1);
phi_hat(n)=phi_hat(n-1)+e(n); % update nCO
end;
% plot waveforms
index_stop=200;
figure;
plot(1:index_stop, phd_output(1:index_stop)),ylabel('Ph. Det.');
figure;
plot(1:index_stop, phi_hat(1:index_stop)*180/pi,'m'),ylabel('Est. Phs.');
index_stop=1000;
% subplot(211)
figure;
plot(1:index_stop, real(nco(1:index_stop)),1:index_stop,real(Signal(1:index_stop))),ylabel('Re-PLL');
% subplot(212)
figure;
plot(1:index_stop, imag(nco(1:index_stop)),1:index_stop,imag(Signal(1:index_stop))),ylabel('Im-PLL');
figure;
plot(1:index_stop,some(1:index_stop));
figure;
plot(1:index_stop,e(1:index_stop));