len=1000;  %仿真數據長度
%QPSK信號源
I_Data=randint(N,1)*2-1;
Q_Data=randint(N,1)*2-1;
s=I_Data + j*Q_Data;
%載波信號
Freq_Sample=2400;   %采樣率
Delta_Freq=-60;     %頻偏
Time_Sample=1/Freq_Sample;   %每個采樣值的持續時間
Delta_Phase=rand(1)*2*pi;   %隨機初始相位
Carrier=exp(j*(Delta_Freq/Freq_Sample*(1:len)+Delta_Phase));  %載波
r=s.*Carrier';  %調制
%鎖相環處理過程如下
Signal_PLL=zeros(len,1);   %鎖相環鎖定及穩定后的數據
NCO_Phase = zeros(len,1);   %鎖定的相位
Discriminator_Out=zeros(len,1);
Freq_Control=zeros(len,1);
PLL_Phase_Part=zeros(len,1);   %鎖相環頻率
PLL_Freq_Part=zeros(len,1);   %鎖相環相位
%環路處理
C1=0.022013;  %環路濾波器參數
C2=0.00024722;
for i=2:N
    Signal_PLL(i)=r(i)*exp(-j*mod(NCO_Phase(i-1),2*pi));   %得到鑒相器的輸入
    I_PLL(i)=real(Signal_PLL(i));   %鑒相器的I路輸入信息數據
    Q_PLL(i)=imag(Signal_PLL(i));   %鑒相器的Q路輸入信息數據
    Discriminator_Out(i)=(sign(I_PLL(i))*Q_PLL(i)-sign(Q_PLL(i))*I_PLL(i))...
                /(sqrt(2)*abs(Signal_PLL(i)));   %鑒相器的輸出
    PLL_Phase_Part(i)=Discriminator_Out(i)*C1;   %環路濾波器處理
    Freq_Control(i)=PLL_Phase_Part(i)+PLL_Freq_Part(i-1);
    PLL_Freq_Part(i)=Discriminator_Out(i)*C2+PLL_Freq_Part(i-1);
    NCO_Phase(i)=NCO_Phase(i-1)+Freq_Control(i);  %進行相位調整
end
%畫圖顯示結果
figure
subplot(2,2,1)
plot(-PLL_Freq_Part(2:len)*Freq_Sample);
grid on;
title('鎖相環頻率響應曲線');
axis([1 len -100 100]);
subplot(2,2,2)
plot(PLL_Phase_Part(2:len)*180/pi);
title('鎖相環相位響應曲線');
axis([1 len -2 2]);
grid on;
%設定顯示范圍
Show_D=300; %起始位置
Show_U=900; %終止位置
Show_Length=Show_U-Show_D;
subplot(2,2,3)
plot(r(Show_D:Show_U),'*');
title('進入鎖相環的數據星座圖');
axis([-2 2 -2 2]);
grid on;
hold on;
subplot(2,2,3)
plot(Signal_PLL(Show_D:Show_U),'r*');
grid on;
subplot(2,2,4)
plot(Signal_PLL(Show_D:Show_U),'r*');
title('鎖相環鎖定及穩定后的數據星座圖');
axis([-2 2 -2 2]);
grid on;
 
figure
%設定顯示范圍
Show_D=600; %起始位置
Show_U=650; %終止位置
Show_Length=Show_U-Show_D;
subplot(2,2,1)
plot(I_Data(Show_D:Show_U));
grid on;
title('I路信息數據');
axis([1 Show_Length -2 2]);
subplot(2,2,2)
plot(Q_Data(Show_D:Show_U));
grid on;
title('Q路信息數據');
axis([1 Show_Length -2 2]);
subplot(2,2,3)
plot(I_PLL(Show_D:Show_U));
grid on;
title('鎖相環輸出I路信息數據');
axis([1 Show_Length -2 2]);
subplot(2,2,4)
plot(Q_PLL(Show_D:Show_U));
grid on;
title('鎖相環輸出Q路信息數據');
axis([1 Show_Length -2 2]);