matlab实现二阶锁相环参数设计

锁相环数学模型

图1 锁相环的复频域模型

锁相环是由鉴相器、环路滤波器、压控振荡器组成的，其中系统传递函数可以表示成$H(s)$，误差传递函数可以表示成$H_e(s)$。其中参考公式可以参考《锁相环技术原理FPGA实现》。在此只给出一些重要的公式。
$$\begin{array}{r}H(s)=\frac{{\theta }_2(s)}{{\theta }_1(s)}=\frac{K_{d}F(s)N(s)}{1+K_{d}F(s)N(s)}\\ H_e(s)=\frac{{\theta }_e(s)}{{\theta }_1(s)}=\frac{1}{1+K_{d}F(s)N(s)}\end{array}$$

二阶锁相环的环路滤波器

环路滤波器是一个低通滤波器，所以我们可以设其$F(s)=k_1+\frac{k_2}{S}$,直接给出$k_1、k_2$的参数设计公式。
$$\begin{array}{r}{k}_1=\frac{8}{3}{B}_n\\ k_2=\frac{1}{2}{k}_1^2\end{array}$$
参考书籍可以参考《Controlled‐root formulation for digital phase‐locked loops》
利用双线性变换公式将$F(s)$转换成相应的Z域表达式$$\begin{array}{r}F(z)=k_1+k_2\frac{T}{2}\frac{1+z^{-1}}{1-z^{-1}}\end{array}$$
所以得到时域的表达形式
$$\begin{array}{r}V(k)-V(k-1)=k_1[ϵ(k)-ϵ(k-1)]+k_{2}T/2[ϵ(k)-ϵ(k-1)]\end{array}$$,$ϵ(k)$是环路滤波器的输入，$V(k)$是环路滤波器的输出。

数字压控振荡器

压控振荡器的Z域表达式在网上到处都可以查到，但是我认为这其实也是最困扰我或者说大部分人的地方。首先还是先给出压控振荡器的Z域模型。
$$\begin{array}{r}N(z)=\frac{K_{0}T{z}^{-1}}{1-z^{-1}}\end{array}$$简单解释一下其中$K_0$是考虑到FPGA实现时候的增益问题，$T$是更新周期这是一个非常重要的参数不同于系统的采样间隔。将$N(z)$转换成时域表达式
$$\begin{array}{r}T\ast V(k-1)={\theta }_0(k)-{\theta }_0(k-1)\end{array}$$如果我想要利用得到的是频率而不是相位则重构公式
$$\begin{array}{r}f(k)=\frac{{\theta }_0(k)-{\theta }_0(k-1)}{2\pi T}\end{array}$$
$$\begin{array}{rl}\delta f(k)& =\frac{{\theta }_0(k)-{\theta }_0(k-1)-{\theta }_0(k-1)+{\theta }_0(k-2)}{2\pi T}\\ & =T\frac{V(k)-V(k-1)}{2\pi T}\\ & =\frac{k_1[ϵ(k)-ϵ(k-1)]+k_{2}T/2[ϵ(k)-ϵ(k-1)]}{2\pi }\end{array}$$到此可以看到一个非常巧妙的地方就是可以直接有环路滤波器的输出直接得到频率的变换率，这也是我这一篇文章想要说明的问题（对于初学者来说可能会对与锁相环怎么去得到频率变化率苦恼，对此给出了公式解答）

matlab代码

// An highlighted block
clc
clear all;

fs = 50e6; %采样频率
ts = 1/fs; 
dataLen = 10e6;  %数据长度
SNR = -15;
realFc = 10000500; %实信号频率
initPhase = 2*pi*realFc*(0:dataLen-1)*ts+pi/4;%输入信号的实际相位，为了对后面得到相位进行比较
data = sin(2*pi*realFc*(0:dataLen-1)*ts+pi/4); %科斯塔斯环的输入信号
pllFc = 10000000; %本地频率

cumulTime = 10000;%累积时间

n = fs/cumulTime; 
nn = [0:n-1];
blockNum = floor(length(data)/n);% 数据块数目
frame = 0;
phase = 0;

T = 1/10000; %锁相环更新时间
cp1 = 8/3*40;   %二阶锁相环参数一，K1，此时等效噪声带宽20Hz
cp2 = cp1*cp1/2*T/2;   %二阶锁相环参数二，K2
vo=0;          %压控振荡器输入
vo_1=0;        %压控振荡器上一个时刻的输入
dfrqFrame(1) = pllFc;
for frame=2:blockNum
expData = exp(1i*(2*pi*pllFc*ts*nn+phase));
sine = imag(expData);   %本地数据
cosine = real(expData);

x = data((1:n)+((frame-1)*n));
 %将数据转换到基带
xSine = x.*sine;
xCosine = x.*cosine;

I = sum(xSine);      %积分累加相当于滤波
Q = sum(xCosine);
pd(frame) = atan(Q/I)/2/pi;   %这里的2pi是为了公式中求delta f做准备

 
%锁相环
deltapd =  pd(frame) -  pd(frame-1);    %解释一下这，例如上一次理论鉴相输出是0.46pi，这一次理论输出是0.92pi，但是
if(abs(deltapd) >= 0.25)               %atan函数的实际输出是-0.08pi，所以差需要修正。类似的和不需要修正因为cp2比较小，cp1较大
     if(deltapd > 0)
        deltapd = deltapd - 0.5;
     else
        deltapd = deltapd + 0.5;
     end
end
%============================环路滤波、delta F求出来=======================
vo = vo_1 + cp1*(deltapd) + cp2*(pd(frame) +  pd(frame-1));  %VCO输入
pllFc = pllFc + (vo - vo_1); %本地复现载波频率更新
vo_1 = vo;
%==========================================================================
dfrqFrame(frame) = pllFc; 
phase = 2*pi*dfrqFrame(frame-1)*ts*n+phase ;   %得到不同块的相位
dphaseFrame(frame) = phase; 

end
figure(1)
plot((0:n:dataLen-1)/fs, realFc*ones(1,length(dfrqFrame)),'r');
hold on
plot((0:n:dataLen-1)/fs,dfrqFrame);
title('锁相环实时频率和输入频率对比')
xlabel('second');
ylabel('freq/hz')
legend('锁相环跟踪','实际的载波频率');

仿真结果分析

经过这种方法设计的锁相环设计参数只有一个环路噪声带宽$Bn$，从而使设计大大的简化。最值得注意的就是PLL锁相环实现频率输出时候的细节，可以参照主要公式。

图2 锁相环实时频率跟踪

图3 锁相环的实时频率跟踪

简单分析一下，图2是压控振荡器更新周期是1.0000e-04，环路噪声带宽是50Hz时候，环路在0.08s左右收敛。图3是压控振荡器更新周期是1.0000e-04，环路噪声带宽是70Hz时候，环路在0.05s左右收敛。环路噪声带宽的选择取决于用户的实际，卫星的载波环大致选择20Hz左右，码环更精细一点，环路噪声带宽大相应的收敛加快，但是带来的捕获跟踪精度会降低。所以设计者需要根据实际进行权衡。