www.gusucode.com > 针对单音信号的科斯塔斯锁相环matlab仿真程序 > 针对单音信号的科斯塔斯锁相环matlab仿真程序/CostasExample1.m

    %% ************************************************************* %%
%                                                                       %
%                 科斯塔斯锁相环示例程序一:单音信号                       %
%                                                                       %
%% ************************************************************* %%
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('锁相环锁定频率情况');


% figure(2),plot(1:LoopN,TmpI,'b',1:LoopN,TmpQ,'r');title('低通滤波后的鉴相信号');
% FFTSig = 20*log10(abs(fftshift(fft(TmpI + 1j*TmpQ))));
% figure(1),plot((-LoopN/2:LoopN/2-1)/LoopN/T,FFTSig);title('接收信号频谱');

% FFTSig = 20*log10(abs(fftshift(fft(PhaseDiff))));
% figure(1),subplot(2,1,1),plot((-LoopN/2:LoopN/2-1)/LoopN/T,FFTSig);title('接收信号频谱');

% temp = 0;
% for i = 1:LoopN
%     VCOIn(i) = c1 * PhaseDiff(i) + temp;    % c1*x(n)+c2*sum(x(1:n-1))
%     temp = temp + c2 * PhaseDiff(i);
% end
% FFTSig = 20*log10(abs(fftshift(fft(VCOIn))));
% figure(1),subplot(2,1,2),plot((-LoopN/2:LoopN/2-1)/LoopN/T,FFTSig);title('接收信号频谱');