suoxinghuan.m

锁相环的仿真以及含有文本

clear all;
clc;close all;
Freq_Sample=1000; %采样率，Hz
Simulation_Time=2; %仿真时间,s
Simulation_Length=Freq_Sample*Simulation_Time; %仿真数据长度
Delta_Freq=10; %频偏，Hz。锁相环中心频率与输入频率的差，这里令锁相环中心频率和初始
%相位为0，只锁定频偏。
Time_Sample=[0:Simulation_Length-1]/Freq_Sample; %采样时间序列
Delta_Phase=rand(1)*2*pi; %随机初相，Rad
Carrier=exp(j*(2*pi*Delta_Freq*Time_Sample+Delta_Phase)); %载波信号
%基带信号
PLL_Mode = 2; %定义锁相环的工作模式：单载波为“1”、BPSK调制为“2”、QPSK调制为“3”
if PLL_Mode == 1
I_Data=ones(1,Simulation_Length);
Q_Data=I_Data;
elseif PLL_Mode == 2
Freq_Bit=50; %bit速率，Hz，采样频率应设为bit速率的整数倍，bit速率最好小于频偏的一半
Point_Number=Freq_Sample/Freq_Bit; %1个bit周期内的采样点数
Bit_Length=ceil(Freq_Bit*Simulation_Time); %准确bit数目可能为小数，取ceil
% Bit_Length_Ori=Freq_Bit*Simulation_Time; %准确bit数目，可能为小数
% Bit_Length=floor(Bit_Length_Ori); %取整
% if Bit_Length_Ori~=Bit_Length
% Bit_Length=Bit_Length+1; %如果有小数，加1
% end
I_Bit=randint(1,Bit_Length)*2-1; %bpsk调制
I_Data=I_Bit(1)*ones(1,Point_Number);
for k=2:Bit_Length
I_Data=[I_Data I_Bit(k)*ones(1,Point_Number)]; %生成各采样点数据
end
Q_Data=zeros(1,length(I_Data));
else
I_Data=randint(1,Simulation_Length)*2-1;
Q_Data=randint(1,Simulation_Length)*2-1;
end
Signal_Source=I_Data + j*Q_Data; %数据bit
%调制处理
Signal=Signal_Source(1:length(Carrier)).*Carrier; %数据的I、Q路都是信号，都通过
%载波发送，不需与载波I、Q相对应

snr=6;
Signal_Channel=awgn(Signal,snr,'measured');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%以下为锁相环处理过程
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%参数清零
Signal_PLL=zeros(1,Simulation_Length);
NCO_Freq = zeros(1,Simulation_Length);
NCO_Phase = zeros(1,Simulation_Length);
NCO_Out=exp(-j*NCO_Phase);
Discriminator_Out=zeros(1,Simulation_Length);
Phase_Control=zeros(1,Simulation_Length);
Freq_Control=zeros(1,Simulation_Length);
PLL_Phase_Part=zeros(1,Simulation_Length);
PLL_Freq_Part=zeros(1,Simulation_Length);
Phase_Diff=zeros(1,Simulation_Length);
%环路参数
gain=2000; %锁相环增益k0*k1
ksi=0.707; %阻尼因子，瞬态响应平坦，不过缓，环路带宽较小
Loopwidth=20; %环路带宽,Hz，与pll稳定性有关
digital_sensitivity=20; %数字灵敏度，k1/fs
k1=digital_sensitivity*Freq_Sample; %nco的增益
k0=gain/k1; %鉴相器后的增益
%以上参数加上采样率即可确定该类锁相环的特性
omega_n=2*Loopwidth/(ksi+1/4/ksi); %自然频率,rad/s
C2=omega_n^2/Freq_Sample/gain; %将有源比例积分滤波器用双线性变换进行数字化即得C1,C2
C1=2*omega_n*ksi/gain-C2/2;

for n=2:Simulation_Length
%鉴相器
Signal_PLL(n)=Signal_Channel(n)*NCO_Out(n-1);
I_PLL(n)=real(Signal_PLL(n));
Q_PLL(n)=imag(Signal_PLL(n));
if PLL_Mode == 1
%四象限反正切
Discriminator_Out(n)=atan2(Q_PLL(n),I_PLL(n))*k0;
% %Q近似
% amplitude_new=abs(Signal_PLL(n));
% if n==2
% amplitude_ave=amplitude_new;
% else
% amplitude_ave=(amplitude_new+amplitude_ave*(n-2))/(n-1);
% end
% Discriminator_Out(n)=Q_PLL(n)/amplitude_ave*k0;
elseif PLL_Mode == 2
%经典Costas鉴别器，加了归一化，使其对幅度不敏感
Discriminator_Out(n)=I_PLL(n)*Q_PLL(n)/abs(Signal_PLL(n))^2*k0;
% %面向判决的Costas鉴别器，加了归一化，使其对幅度不敏感
% sgn=cal_sign(I_PLL(n));
% Discriminator_Out(n)=sgn*Q_PLL(n)/abs(Signal_PLL(n))*k0;
% %次最佳Costas鉴别器
% Discriminator_Out(n)=Q_PLL(n)/I_PLL(n)*k0;
% %二象限反正切
% Discriminator_Out(n)=atan(Q_PLL(n)/I_PLL(n))*k0;
else
sgn=cal_sign(I_PLL(n));
Discriminator_Out(n)=(sgn*Q_PLL(n)-sgn*I_PLL(n))/(sqrt(2)*abs(Signal_PLL(n)))*k0;
end
%有源比例积分滤波器
PLL_Phase_Part(n)=Discriminator_Out(n)*C1; %暂态相差
%Discriminator_Out是相位差 
PLL_Freq_Part(n)=Discriminator_Out(n)*C2+PLL_Freq_Part(n-1); %稳态相差，C2可视为暂态相差对稳态相差的加权系数
Phase_Control(n)=PLL_Phase_Part(n)+PLL_Freq_Part(n); %经过滤波的一个采样间隔内的相位差

%数控振荡器，中心频率和初始相位设为0。实际只能调频，通过调频来调相
% NCO_Phase(n)=NCO_Phase(n-1)+Freq_Control(n)/Freq_Sample*k1;
% NCO_Out(n)=exp(-j*NCO_Phase(n));
Phase_Diff(n)=Phase_Control(n)/Freq_Sample*k1; %加权后一个采样间隔的相位差
NCO_Freq(n)=Phase_Diff(n)*Freq_Sample; %除以采样间隔，得到pll中心频率与输入信号的频率差
NCO_Phase(n)=NCO_Phase(n-1)+NCO_Freq(n)/Freq_Sample; %本地载波相位
NCO_Out(n)=exp(-j*NCO_Phase(n)); %本振信号
end
PLL_Freq=PLL_Freq_Part*Freq_Sample/2/pi*digital_sensitivity; 
%Hz,锁定的频率为归一化角频率，即freq_pll=f0/fs*2*pi
PLL_Phase=Discriminator_Out*180/pi; %度 %PLL_Phase_Part(2:Simulation_Length)*180/pi;
%画图显示结果
figure
subplot(2,2,1),plot(Time_Sample,PLL_Freq,'*')
grid on
xlabel('时间(s)')
ylabel('输出频率(Hz)')
title('锁相环频率响应曲线')
%axis([1 Simulation_Length -100 100])
subplot(2,2,2),plot(Time_Sample,PLL_Phase)
grid on
xlabel('时间(s)')
ylabel('相位差(°)')
title('锁相环相位响应曲线')
%axis([1 Simulation_Length -2 2])

%设定显示范围
Show_D=1000; %起始点序号
Show_U=1800; %终止点序号
Show_Length=Show_U-Show_D;
subplot(2,2,3),plot(Signal_Channel(Show_D:Show_U),'*') %信号点星座图
%此处plot函数画的每个点以其实部为横坐标，虚部为纵坐标
grid on
title('进入锁相环的数据星座图')
axis([-2 2 -2 2])
hold on %使得下面对于子图3的处理与上面的叠加起来
subplot(2,2,3),plot(Signal_PLL(Show_D:Show_U),'r*') %载波剥离后的输出信号
%锁相环的工作过程应该是先用NCO输出的本地信号下变频，然后鉴相、滤波、调整NCO。当处于锁定
%状态，下变频后只有数据bit信号和噪声，可以视为锁相环的输出，但更应该视为用锁相环锁定的频率
%进行载波剥离后的输出信号
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=1000; %起始点序号
Show_U=1800; %终止点序号
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]);