close all;
clear all;
clc;
%% 信号及锁相环参数
Fs = 1e6; % 采样率
F0 = 1e5; % 信号载波
DetaF = 100; % 本地时钟与信号载波的差,频差越大需要的锁相周期越小
Time = 1; % 信号总时长
snr = 10; % 信噪比
T = 1e-3; % 锁相周期,要求远小于1/DetaF
LoopN = Time / T; % 锁相次数
Len = Fs*T; % 每次锁相采用的采样点数
SigLen = Fs*Time; % 信号总长
wfc = 2*pi*(F0+DetaF); % 初始本振角频率
t = 0:1/Fs:T-1/Fs; % 各段锁相中的时间序列
RecSig = cos(2*pi*F0*(0:SigLen-1)/Fs); % 接收端信号
RecSig = awgn(RecSig,snr);
% FFTSig = 20*log10(abs(fftshift(fft(RecSig))));
% figure(1),subplot(2,1,1);plot((-SigLen/2:SigLen/2-1)/SigLen*Fs,FFTSig);title('接收信号频谱');
% figure(1),subplot(2,1,2);plot(RecSig);title('接收信号时域波形');
%% 锁相:鉴相+环路滤波+VCO
c1 = 153.7130; % 环路滤波器的参数
c2 = 6.1498;
phase = 0;
temp = 0;
% 记录中间过程的数组
PhaseDiff = zeros(1, LoopN);
FrqDiff = zeros(1, LoopN);
TmpI = zeros(1, LoopN);
TmpQ = zeros(1, LoopN);
VCOIn = zeros(1, LoopN);
for i = 1:LoopN
% 产生本振参考信号:根据当前锁相后的角频率
LoOsc = exp(1j*(wfc*t+phase));
sine = imag(LoOsc);
cosine = real(LoOsc);
% 鉴相:输出本振与接收信号相位差
% 1、接收信号与本地参考信号相乘
x_sine = RecSig((i-1)*Len+1:i*Len) .* sine;
x_cosine = RecSig((i-1)*Len+1:i*Len) .* cosine;
% 2、低通滤波+抽取(降低处理速率),滤波器为全1的门函数,对应频域为sa函数
TmpI(i) = sum(x_cosine);
TmpQ(i) = sum(x_sine);
PhaseDiff(i) = atan2(TmpQ(i), TmpI(i)); % 得到锁相环的输入
% 环路滤波器:产生VCO输入
VCOIn(i) = c1 * PhaseDiff(i) + temp; % c1*x(n)+c2*sum(x(1:n-1))
temp = temp + c2 * PhaseDiff(i);
% VCO
wfc = wfc - VCOIn(i) * 2 * pi; % 改变本地频率
FrqDiff(i) = wfc; % 记录频率变化
phase = wfc * Len / Fs + phase; % 得到不同块的相位
end
figure(3),plot(1:LoopN,FrqDiff/(2*pi),'b',1:LoopN,F0,'r');
legend('锁相环跟踪频率','实际的载波频率');title('锁相环锁定频率情况');
23