第4章（2）内容如下:
- 一、用二阶锁相环对有多普勒频偏和载波随机相位的psk信号进行相干解调
- 二、总结
上一讲《载波同步与数字锁相环仿真（1）》主要讲解了一下锁相环的基本原理及其MATLAB实现，观察锁相环是如何跟踪信号相位，并且最终可以和接收到的信号同频同相。
其中，二阶锁相环的核心代码需要你扎实学习与掌握。
理解了锁相环的基本原理后，仍然有一些问题需要我们去解决，比如：
1、如何用数字锁相环对有多普勒频偏的2psk信号进行相干解调？
2、锁相环进入同步状态的时间要多久呢？能否调整呢？锁相环会不会跟不上信号的相位变化呢？
3、如果接收到的信号包含噪声，锁相环能否正常工作呢？或者说锁相环正常工作条件是什么呢？
一、用二阶锁相环对有多普勒频偏和载波随机相位的psk信号进行相干解调
我刚开始好不容易搞懂了二阶锁相环的MATLAB实现，写在了《数字锁相环的MATLAB实现》，而如何将其运用到有多普勒频偏的2psk信号的相干解调，调试代码与整理思路又花了好些天时间。
在写有多普勒频偏的2psk信号的相干解调时，我便充分认识到自己对于锁相环理论知识的欠缺。检查很多遍代码后，发现代码框架没有大问题，但就是死活不出来正确结果，这样的心情真是很磨人。
于是我从图书馆借出来西电郑继禹的《锁相技术》和杜勇的《锁相环原理及FPGA实现》，进行学习。然而理论知识太多，我还未融会贯通。
回到有多普勒频偏的2psk信号的相干解调，上一讲《载波同步与数字锁相环仿真（1）》的代码中，进入锁相环的信号还没有乘以双极性序列，即没有对基带信号进行调制。
我在《BPSK调制解调器仿真》的第三点“如何保证仿真正确性”中讲过：测试每个模块的正确后，再加入到大程序仿真中。
即当我们写一个大的程序时，为了避免总程序出错不知道从什么地方开始找，我们在最开始写每个程序模块的时候，先对其进行测试，测试正确后，加入到大程序后，来减少程序调试的错误。
这次将双极性序列乘以有多普勒频偏的载波进行测试。
注意到，由于发送端乘以 cos w t \cos wt coswt ,接收端乘以 cos w t \cos wt coswt +低通滤波，会把2w的频率滤除，只剩下有多普勒频偏的基带频率。因此，直接乘 cos ( 2 p i ⋅ f d o p ⋅ t + p h i 0 ) \cos \left( {2pi \cdot {f_{dop}} \cdot t + phi0} \right) cos(2pi⋅fdop⋅t+phi0) ，验证程序正确后，之后再乘以 cos ( 2 p i ⋅ ( f d o p + f c ) ⋅ t + p h i 0 ) \cos \left( {2pi \cdot \left( {{f_{dop}} + {f_c}} \right) \cdot t + phi0} \right) cos(2pi⋅(fdop+fc)⋅t+phi0) ，后面的程序将看到这个过程。
%%%%%%%%%%%%%%%%%%%%% 环路滤波器原理仿真文件 %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%% second_loopfilter_sim4.m %%%%%%%%%%%%%
%%%%%%%%%%%%%%%% data:2020年10月2日 author:飞蓬大将军 %%%%%%%%%%
%%%%%程序说明
%%%利用锁相环完成频率跟踪过程的原理仿真，采用一个频率固定的信号源
%无加噪声过程
%
%%% sim系列说明
%sim3：未乘以1，-1序列
%sim4：%乘以1，-1序列
%%%%%% 仿真环境
%软件版本：MATLAB R2019a
clc;
clear all;
close all;
fdop = 30; %多普勒偏移50Hz
phi0 = pi/6;
fs = 1e4; %采样频率
T = 1/fs;
n = 3000;
tend = n*T;
t1 = 0:T:tend - T; %1到3000
phase1= 2*pi*fdop.*t1 + pi/6; %输入信号相位变化
phase2 = zeros(1,length(t1));
frame_length = 20;
frame = ones(1,frame_length);
data =[frame randi([0,1],1,n-frame_length)];
bipolar_data = 2*data-1;
I_signal = bipolar_data.*cos(2*pi*fdop.*t1 + phi0);
Q_signal = bipolar_data.*sin(2*pi*fdop.*t1 + phi0);
%%%%加入噪声
signal = I_signal + 1j*Q_signal;
signal_temp = awgn(signal,30,'measured');
%%%%%%%%%%%%环路滤波器参数设置
%%%%%2阶锁相环参数
Bn = 400; %噪声等效带宽
ts = 1/fs;%时间周期
damp = 0.707;%衰减系数
Kp = 1; %K0*K1相乘，下面公式默认了k0 = k1 =1
Wn = 2*Bn/(damp+1*(4*damp));
plus = 1; %环路滤波器增益
C1 = plus*8*damp*Wn*ts/(Kp*(4+4*damp*Wn*ts+(Wn*ts)^2));%环路参数1，C1 = alpha
C2 = plus*4*(Wn*ts)^2/(Kp*(4+4*damp*Wn*ts+(Wn*ts)^2)); %环路参数2 C2 = beta
%%%%%采用等效变换积分器的2阶锁相环
A = zeros(1,n);
B = zeros(1,n);
error = zeros(1,n);
for k = 1:n
rotate_sig = signal_temp(k)*exp(-1j*phase2(k));
%%%分I、Q两路
sig_I_pll(k) = real(rotate_sig);
sig_Q_pll(k) = imag(rotate_sig);
%%%-------松尾环鉴相算法
u1(k) = sign(sig_I_pll(k));
u2(k) = sign(sig_Q_pll(k));
e(k) = sign(sig_I_pll(k))*sig_Q_pll(k)/sqrt(sig_I_pll(k)^2 + sig_Q_pll(k)^2);
if (k>=2)
%%%环路滤波器
A(k) = A(k-1) + C2*e(k) ;
B(k) = A(k) + C1*e(k);
%%%相位累加
phase2(k+1) = phase2(k)+ B(k);
end
%%%%%%为了看到频差
if(k>200) &&(k<n-200)
freq(k-200) = (phase2(k)-phase2(k-100))/(2*pi)/100/ts;
end
end
u1_temp = (u1+1)/2;
[err_number,bit_err_ratio]= biterr(data,u1_temp);
x = (1:n)/fs;
figure(1);
plot(x,e);
title('二阶锁相环，初始频差500Hz，噪声带宽300Hz');
xlabel('时间（s）');
ylabel('相差');
grid on ;
figure(2);
plot(freq);
title('二阶锁相环，初始频差50Hz，噪声带宽300Hz');
xlabel('时间（s）');
ylabel('频差');
%%%%%%%%%%%%%%% 实验结论
%%%%能够跟踪信号的频率
图1 相差随时间变化图
图2 频差随时间变化图
在上面的代码中，有极性鉴相器和松尾环鉴相器，这部分内容可以参考知网论文《卫星定时接收机的关键技术研究》，里面对捕获和跟踪有介绍，值得仔细阅读。
在《载波同步与数字锁相环仿真（1）》和《数字锁相环的MATLAB实现》的代码中，我们得出结论：噪声带宽值增加，捕获速度会加快，但是抖动会增加。
为了在一个符号内锁相环能迅速捕获并跟踪到信号，后续的抖动便会增加。对应到上面second_loopfilter_sim4.m 文件，读者可调整Bn也可进行验证。
有了多普勒频偏，怎么来设置相应的Bn？有没有什么规律呢？这个问题我还在思考中。
接下来，说一下噪声对锁相环路的影响。
详细内容可以看张厥盛的《锁相技术》第三章，本文不加证明的直接摘录书中一些重要结论，方便后续的代码实现。
锁相环的噪声与干扰主要是两个来源：一类是与信号一起进入环路的输入噪声与谐波干扰；另一类是环路部件产生的内部噪声与谐波干扰，以及压控振荡器控制端感应的寄生干扰等，其中压控振荡器内部噪声是主要的噪声源。
想一想，说这两类噪声与干扰有什么用呢？
在解决问题时，为了达到所有的优化目标往往是不可能的，常常需要trade-off。减小环路噪声等效带宽，可以有效滤除第一类噪声，但是压控振荡器的相位噪声主要集中在低频部分，所以这会使得压控振荡器的相位噪声大量通过。
就我个人认为，因为在仿真的过程中，我们不涉及压控振荡器的相位噪声仿真，所考虑的噪声只有第一类，也就是对信号加上的高斯白噪声。
如果这样的话，我们是不是可以使得环路噪声等效带宽设置尽可能小呢？
当然不行，因为等效噪声带宽又会影响捕获带和抖动情况。
有了上面的这段代码，我们怎么知道解调是否正确呢？
与理论的误码率曲线进行对比，是不错的方法，因此将实验所得的误码率与berawgn函数进行对比。
%%%%%%%%%%%%%%%%%%%%% 环路滤波器原理仿真文件 %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%% second_loopfilter_sim6.m %%%%%%%%%%%%%
%%%%%%%%%%%%%%%% data:2020年10月2日 author:飞蓬大将军 %%%%%%%%%%
%%%%%程序说明
%%%利用锁相环完成频率跟踪过程的原理仿真，采用一个频率固定的信号源
%%% sim系列说明
%sim3：未乘以1，-1序列
%sim4：乘以1，-1序列
%%%%%% 仿真环境
%软件版本：MATLAB R2019a
clc;
clear all;
close all;
%%%%%%%%%%%%%%%%%%%%%%% 程序主体 **************************
fdop = 20; %多普勒偏移20Hz
phi0 = pi/6;
fs = 1e4; %采样频率
T = 1/fs;
% n = 3000;
n = 10000000;
tend = n*T;
t1 = 0:T:tend - T; %1到3000
phase1= 2*pi*fdop.*t1 + pi/6; %输入信号相位变化
phase2 = zeros(1,length(t1));
frame_length = 20;
frame = ones(1,frame_length);
data =[frame randi([0,1],1,n-frame_length)];
bipolar_data = 2*data-1;
% I_signal = bipolar_data.*cos(2*pi*fdop.*t1 + phi0);
% Q_signal = bipolar_data.*sin(2*pi*fdop.*t1 + phi0);
%
%
% %%%%加入噪声
% signal = I_signal + 1j*Q_signal;
signal = bipolar_data.*exp(1j*(2*pi*fdop.*t1 + phi0));
signal_temp = awgn(signal,9,'measured');
%%%%%%%%%%%%环路滤波器参数设置
%%%%%2阶锁相环参数
Bn = 200; %噪声等效带宽
ts = 1/fs;%时间周期
damp = 0.707;%衰减系数
Kp = 1; %K0*K1相乘，下面公式默认了k0 = k1 =1
Wn = 2*Bn/(damp+1*(4*damp));
plus = 1; %环路滤波器增益
C1 = plus*8*damp*Wn*ts/(Kp*(4+4*damp*Wn*ts+(Wn*ts)^2));%环路参数1，C1 = alpha
C2 = plus*4*(Wn*ts)^2/(Kp*(4+4*damp*Wn*ts+(Wn*ts)^2)); %环路参数2 C2 = beta
%%%%%采用等效变换积分器的2阶锁相环
A = zeros(1,n);
B = zeros(1,n);
error = zeros(1,n);
for k = 1:n
rotate_sig = signal_temp(k)*exp(-1j*phase2(k));
%%%分I、Q两路
sig_I_pll(k) = real(rotate_sig);
sig_Q_pll(k) = imag(rotate_sig);
%%%-------松尾环鉴相算法
u1(k) = sign(sig_I_pll(k));
u2(k) = sign(sig_Q_pll(k));
% u3(iii) = sign(sig_I_pll(iii)+sig_Q_pll(iii));
% u4(iii) = sign(sig_I_pll(iii)-sig_Q_pll(iii));
% e(k) = sign(u1(iii)*u2(iii)*u3(iii)*u4(iii));
%%%%%BPSK极性鉴相算法
e(k) = sign(sig_I_pll(k))*sig_Q_pll(k)/sqrt(sig_I_pll(k)^2 + sig_Q_pll(k)^2);
if (k>=2)
%%%环路滤波器
A(k) = A(k-1) + C2*e(k) ;
B(k) = A(k) + C1*e(k);
%%%相位累加
phase2(k+1) = phase2(k)+ B(k);
end
%%%%%%为了看到频差
if(k>200) &&(k<n-200)
freq(k-200) = (phase2(k)-phase2(k-100))/(2*pi)/100/ts;
end
end
u1_temp = (u1+1)/2;
[err_number,bit_err_ratio]= biterr(data,u1_temp);
x = (1:n)/fs;
figure(1);
plot(x,e);
title('二阶锁相环，初始频差20Hz，噪声带宽300Hz');
xlabel('时间（s）');
ylabel('相差');
grid on ;
figure(2);
plot(freq);
title('二阶锁相环，初始频差20Hz，噪声带宽300Hz');
xlabel('时间（s）');
ylabel('频差');
%%%%%%%%%%%%%%%%%%%%%%%% 理论值 %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% EbN0(dB) berawgn(EbN0,'psk',2,'nondiff')
%%%% 3 0.022878407561085
%%%% 4 0.012500818040738
%%%% 5 0.005953867147779
%%%% 6 0.002388290780933
%%%% 7 0.000772674815378
%%%% 8 0.000190907774076
%%%% 9 0.000033627228420
%%%% 10 0.000003872108216
%%%%%%%%%%%%%%%%%%%%%%%% 实验结果 %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% EbN0(dB) bit_err_ratio
%%%% 6 0.002394300000000
%%%% 7 7.739999999999999e-04
%%%% 8 1.930000000000000e-04
%%%% 9 3.390000000000000e-05
%%%% 10 3.900000000000000e-06
%%%%%%%%%%%%%%% 实验结论
%%%%能够跟踪信号的频率并正确解调出数据
通过上面的代码，已经验证过锁相环可以用于有多普勒频偏的2psk信号的相干解调。读者也可继续测试不同的 E b / N 0 {E_b}/{N_0} Eb/N0 下，得到的误码率。
上面的代码，没有基带成型、匹配滤波、上载波、捕获、下载波、信道编解码等过程。结合《第1章：BPSK调制解调器仿真》和《第2章：线性分组码》，便可以把代码补充完整。
至此，用锁相环进行载波同步的过程已经介绍完毕。
那么多普勒频偏达到多少，锁相环便会跟不上呢？即锁相环无法进入锁定状态。
举个例子来说，假如多普勒频偏10Hz，锁相环可以快速跟上，但是几百Hz频偏，锁相环便不能快速跟上了，需要的时间可能会很长。
在卫星通信中，由于卫星相对地面终端的移动速度特别快，多普勒频偏便会很大。
那有没有办法，让进入锁相环的信号频偏比较小呢?
在进入锁相环之前，可对信号进行初步的频偏补偿。举个例子，假如频偏是220Hz，我们用0-100Hz，100-200Hz，200-300Hz这样一段一段的频率去扫描。然后惊奇发现，频偏范围是属于200-300区间的，于是先补偿200Hz给信号，再进入锁相环，对信号进行跟踪。
那怎么判断出来频偏范围是200-300Hz这个区间的呢？
就是看信号频率离谁越近，就认为是落在哪个区间。因此，本质是最大似然估计的思想。
以上的用0-100Hz，100-200Hz，200-300Hz这样一段一段的频率去扫描的过程，其实信号捕获中频偏粗估计的过程。
值得注意的是，信号捕获和锁相环里面的捕获不是一个概念，信号捕获是指接收端判断信号有没有来了，来了的话，就是捕获成功，没来的话，就是捕获失败。
二、总结
从我个人调试代码的过程中发现，锁相环虽然好用，但是参数设置比较麻烦，理论分析又需要非常非常非常扎实的理论功底，而我目前还没有这么强的扎实功底。
因此，我是结合上述书籍与材料中的理论分析，然后不断的调整参数，比如环路噪声等效带宽Bn，环路增益等，来使得锁相环满足跟踪频率和解调信号的要求。这也再一次体现“理论指导实践，实践印证理论”的过程。
锁相环的内容很丰富，但是本科一般都不学习，确实比较难。在自学这部分的内容时，我对阅读又有了两点体会：
1、读经典材料
我现在越来越明白为什么大家说要读经典书籍和论文了，经得起时间考验的东西就是不太一样。比如我在上面提到的张厥盛的《锁相技术》，这可能会比你随便找一本锁相理论的书阅读起来舒服。
2、相互印证
读书籍或者论文，也别只读一本。最好是相同的知识点，看不同的书籍怎么讲，谁讲的对，谁讲的好，高下立判。虽然多花了一些时间，好像在阅读相同的材料，但是不同的人讲解知识的角度、侧重点，这也会让我有所收获，学习便也扎实了些。
特别感谢微信公众号“通信工程师专辑”的蔡凡博士，本人虽学习兴趣浓厚，但能力有限，且缺少大型工程项目实践，不足甚至错误之处，真的在所难免，欢迎你批评指正。
也欢迎读者朋友就以上讲到的相关技术问题与我交流，一起学习，共同进步。请你也别忘了把这篇文章分享给你身边正在学习通信专业的同学们，也许能够帮到Ta。这是《陈老湿·通信MATLAB》仿真的第4章，期待下次更新见！