% 本程序仿真正弦波的跟蹤過程
% PLL illustration using MATLAB 
% 2002
clear all；
% close all；
% 定義初始相位偏移和輸入連續波頻率以及采樣頻率
theta = 60*pi/180; 
f=1e3; 
fs=100e3;

% 生成未調制連續波的實部和虛部
k=1:1:1000;
delf=f/20;
cpx1=exp(j*(2*pi*k*(f+delf)/fs+theta))+.01*(rand(1,1000)+j*rand(1,1000));


%初始化鎖相環 
phi_hat(1)=30;
e(1)=0;
phd_output(1)=0;
nco(1)=0

% 定義環路濾波器參數
kp=0.15;        % 比例常數
ki=0.1;            % 積分常數

% 鎖相環的實現
for n=2:length(cpx1)
    nco(n)=conj(exp(j*(2*pi*n*f/fs+phi_hat(n-1))));                     % 數控振蕩器
    phd_output(n)=imag(cpx1(n)*nco(n));                                  % 鑒相
    e(n)=e(n-1)+(kp+ki)*phd_output(n)-ki*phd_output(n-1);     %濾波
    some(n)=(kp+ki)*phd_output(n)-ki*phd_output(n-1);
    phi_hat(n)=phi_hat(n-1)+e(n);                                               % 更新數控振蕩器
end;

% 繪圖
index_stop=200;
figure
subplot(211),plot(1:index_stop, phd_output(1:index_stop)),ylabel('Ph. Det.')
subplot(212),plot(1:index_stop, phi_hat(1:index_stop)*180/pi),ylabel('Est. Phs.')
figure, index_stop=200;
subplot(211),plot(1:index_stop, real(nco(1:index_stop)),1:index_stop, real(cpx1(1:index_stop))),
ylabel('RE-PLL')
subplot(212),plot(1:index_stop, imag(nco(1:index_stop)),1:index_stop, imag(cpx1(1:index_stop))),
ylabel('IM-PLL')