数字锁相环MATLAB代码
奈奎斯特型全数字锁相环（NR-DPLL）
注：本文截取于通信原理课程综合设计，载波提取部分中的锁相环解调部分中的基础锁相环。MATLAB编程仿真实现，想要simulink实现的同学要失望啦。代码在本文末，抱歉未加注释。理解本文需要的知识：信号与系统，数字信号处理，同步技术。
2.7载波的同步提取
提取载波信息可用锁相环进行跟踪载波或调制信息。本文采用奈奎斯特型全数字锁相环（NR-DPLL ）对接收信号进行载波同步提取，并用于相干解调。 2.7.1 NR-DPLL 结构介绍
数字锁相环的基本组成如下
图2-6 数字锁相环的组成
NR-DPLL 是基于奈奎斯特采样鉴相器、数字环路滤波器、数字控制振荡器的一种数字锁相环。下面分别对各部分作简要介绍。 2.7.2 奈奎斯特采样鉴相器
奈奎斯特采样鉴相器的组成框图如图2-7所示。
图2-7 奈奎斯特采样鉴相器的组成框图
为了表述方便，设数字控制振荡器（NCO ）输出的本振数字信号为
0002()cos(())k k k u t U t t ωθ=+ (2.7-1) 输入信号
101()sin(())i u t U t t ωθ=+ (2.7-2)
其中
100()(),i i o t t t θωθωωω=?+?=-
输入信号经A/D 采样后，第k 个采样时刻采样量化后的数字信号为
01()sin(())i k i k k u t U t t ωθ=+ (2.7-3) 对输入信号进行A/D 变换的采样速率由带通信号奈奎斯特采样定理确定，但为防止信号频谱混叠并保证信号相位信息的有效抽取，采样速率一般选取前置带通滤波器的两倍带宽以上。
令()(),()()i k i o k o u t u k u t u k ==，即()i u k 和()o u k 相乘后，经低通滤波得到的数字误差信号
()sin ()d d e u k U k θ= (2.7-4)
式中
12()()()e k k k θθθ=- (2.7-5)
2.7.3 数字环路滤波器
数字环路滤波器与模拟环路中环路滤波器的作用是一样的，都是为了抑制高频分量及噪声，且滤波器的参数直接影响环路的性能。在实际应用中一阶数字环路滤波器的实现形式如图2-8所示。
图2-8 一阶数字环路滤波器的实现形式
其Z 域传递函数：
2
11
()z ()1c d u k G F G u k z
-=+-（）=
(2.7-6) 按照图2-8中所实现的数字滤波器，其频率特性与理想积分滤波器的频率特性一致；两种滤波器参数之间也有着一定的对应关系。
对理想积分滤波器的传递
函数式采用双线性变换，即令1
1
211s z s T z ---=+，得到
22
1
11
1121()211s s T T G F z G z z τττ---=
+=+-- (2.7-7) 式中,s T 为采样周期，
21122s
T G ττ-= (2.7-8) 21
s
T G τ=
(2.7-9)
式(2.7-8)和式(2.7-9)表明了两种滤波器参数之间的对应关系，也可以说明图2-8所示的一阶数字滤波器就是模拟理想积分滤波器的数字化表示形式。
二阶数字滤波器可由两个一阶数字滤波器串联得到。 2.7.4 数字控制振荡器（NCO ）
NCO 采用直接数字相位综合技术（DDS ），该技术主要是由时钟驱动读取三角函数表，功能框图如图2-8所示。
图2-9基于DDS 的NCO 结构
Ok,截取结束。接下来是编程中必须的参数计算以及思路分析。 首先是参数选择，时域响应最好取
ξ=0.707，
根据经验
12=10ττ，
令NCO 相位控制增益
0K =1
’，鉴相灵敏度d K 1=，则0K s f =，环路增益
d 0K s
K K f ==,
再由理想二阶环路中
2
2
n
τξω=
，
n ω=
于是，
21
1
200
S
T G τ=
=
, 221120.10.097522
S T G
G ττ-=
=-= 以上两个参数同时缩小相同的倍数，捕获带也缩小相同的倍数。而扩
大时，最大到两倍，大于两倍时捕获带不再变化，只是捕获时间短一些。
环路的3Db 带宽
2c n ωΩ==接下来是编程思想，离散信号过系统实际是系统冲击响应与信号卷积，所谓卷积就是错位相加，于是引申出一种叫做重叠保留法的计算方法，即信号可分段计算。鉴相器输入的两个信号每次分别都仅是一个点，相乘后进入滤波器，此处就得用到这个所谓的重叠保留法。输入一个点输出一个点，其他点保留与下次错位相加，再输出一个点。
鉴相器输出的信号过环路滤波器，通过Z 域传递函数计算。 最后是DDS ，网上有很多相关文档，自己看看，注意输入相位要累积，所谓驱动嘛，当然每次依然是输出一个点，这个点作为下次鉴相器输入信号之一。
仿真结果如下：
上图中输入频率是27Hz，锁相环自然频率为32Hz。
搞懂了这些编程思想，以及DDS后，自然谐振频率，输入信号频率自己去做调整。但注意锁相环的参数限制。
到这里，NR-DPLL就完成了。在此基础上平方环，科斯塔斯环，判决反馈环等等都很容易编程仿真实现。
仿真代码：
锁相环：
% function uo=pll(ui,Fs)
clear all
clc
fo=32;
fi=27;
% fi=length(ui)*1/(n*fo);
g1=0.005*2;
g2=0.0975*2;
n=64;
t=0:1/(n*fo):0.5-1/(n*fo);
ui=sin(2*pi*fi*t);%输入信号
uo=zeros(1,length(ui));
uo(1)=dds_sin(10,n*fo,fo,0,1);%Nco输出初始信号
uo1(1)=dds_sin(10,n*fo,fo,0,1);
% erro=ph(ui(1),uo(1));%鉴相器输出
sig(1)=ui(1)*uo(1);
l=30;
hn=fir1(l,0.2*fo/(n*fo));
temp=conv(hn,sig(1));
temp1=temp(2:end);
erro(1)=temp(16);
lpu1(1)=erro(1)*(g1+g2);
lpu2(1)=lpu1(1)*(g1+g2);
m=2*pi*fo/2048*(1:length(t));
for i=1:length(ui)-1
uo(i+1)=dds_sin(10,n*fo,fo,m(i)+lpu2(i)*1,1);
uo1(i+1)=dds_sin(10,n*fo,fo,m(i)+lpu2(i)*1-m(14),1);
sig(i+1)=ui(i+1)*uo(i+1);
temp=conv(sig(i+1),hn);
temp(1:l)=temp(1:end-1)+temp1;
temp1=temp(2:end);
erro(i+1)=temp(16);
lpu1(i+1)=g1*erro(i+1)+g2*erro(i+1)+g2*erro(i);
lpu2(i+1)=g1*lpu1(i+1)+g2*lpu1(i+1)+g2*lpu1(i)+lpu2(i);
end
% end
DDS代码
function [A] = dds_sin(N,fs,f,phi0,Nig)
% % [返回参数] (传递参数)
% f=8;
% N=8; %频率信号字长N:表示频率信号样点的字长
% fs=256; % 采样频率fs:DDS内部采样频率(输出频率信号的采样频率)
% phi0=0; % 初始相位phi0:输出信号的初始相位
% % df=1; %最小频率间隔df:DDS输出信号的最小频率间隔(频谱分辨率)
% Nsig=256; %输出频率信号样点个数Nsig
% % f0 输出信号频率
Nperiod=2^N; %这么多样点表示一个周期正弦信号
n=0:(Nperiod-1);
Sroot=sin(2*pi*n/Nperiod);%*(2^(N-1)-1); Asin=Sroot;
dphi=2*pi/Nperiod;
M=floor(2^N*f/fs);
M0=floor(phi0/dphi);
% A=zeros(1,Nsig);
% A=Asin(x);
for i=1:Nig
x=mod(M0+1+(i-1)*M,Nperiod);
if x==0
x=1;
end
A(i)=Asin(x);
end
% A1=A;
end
% x=mod(M0+(i-1)*M,Nperiod);
% xA=0:1/fs:(Nsig-1)/fs;
% plot(xA,A)
% title('DDS正弦信号时域图');xlabel('t/s'); % end