www.pudn.com > spll_simplest_IQ.rar > spll_simplest_IQ.m, change:2011-08-01,size:5541b


%% 软件锁相环(科斯塔斯环最简单模型) 
 
%双线性变换法% 
%乘法鉴相器% 
%双路科斯塔斯环% 
 
%%  
clear all; 
% close all; 
format long; 
%% -----产生输入正弦信号------ 
 
F0 =1e3;                         %正弦信号频率 
Fs =10e3;                        %采样频率 
  %计算VCO的固有频率 
frag1=fix(F0/(Fs/2)); 
frag2=rem(F0,(Fs/2)); 
if rem(frag1,2)==0&&frag2>=0  %平移 
    Fn=F0-(Fs/2)*frag1; 
else 
    Fn=(Fs/2)*(frag1+1)-F0;%反折 
end 
     
num=10000;%采样点数 
Ts=1/Fs; 
T =num*Ts;%采样的时间长度,单位:秒 
 
%固定多普勒频偏 
Fpian=100;  
 
t = Ts:1/Fs:T;             %产生时间序列 
P0 =0;     
Signal=sin(2*pi*(F0+Fpian).*t+P0); 
 
%---加入噪声---- 
snri=30; 
snr = 10^(snri/10); 
segma = sqrt(1/(2*snr)); 
noise = normrnd(0,segma,1,length(Signal));%产生正态分布噪声 
Signal = Signal+noise;   %加入噪声 
% showfft(Signal,Fs,'输入信号频谱'); 
 
%% ---环路参数设计----- 
 
S=sqrt(2)/2;  %阻尼系数 
 
% BL=(1+4*S^2)/(16*S^2)*max(Fpian);  %     1+4*S^2                 3                                              % 
%                                    % BL=---------- * max(Fpian)=--- * max(Fpian) 保证最大多普勒频偏在快捕带以内 %  
%                                    %      16*S^2                 8                                              %  
 
BL=20; 
 
fn=8*S/(1+4*S^2)*BL; 
wn=2*pi*fn; 
%通过定值BL求快捕带 
dw_able=2*fn*S;  %单位:Hz 
 
%双线性变换法的系数,只受三个参数(wn,S,Ts)的影响 
K1=2*S/wn+Ts/2; 
K2=2*S/wn-Ts/2; 
K3=wn^2*Ts/2; 
A=1;% 环路滤波器的修正系数 
 
%%   %%%%%%%%%%%%理论性能指标%%%%%%%%%%%%%%% 
 
BL_ini=wn*(1+4*(S^2))/(16*pi*S);%环路等效噪声带宽,单位:Hz 
dw_ini=2*wn*S/(2*pi); %环路快捕带 
 
Fpian_max=max(Fpian);%最大多普勒频偏 
 
df0=Fpian_max; %环路输入信号与本振的初始频差 
dw0=2*pi*df0; 
Tp=dw0^2/(2*S*(wn^3));%环路捕获时间,单位:s 
ts=4/(S*wn);%调节时间,单位:s 
SNRL=8.62+10*log10(BL_ini);%环路锁定信噪比门限,单位:dB 
   Bi=0;%环路输入信号的有效带宽  
SNR0=snri+10*log10(Bi/2*BL_ini);%环路输出信噪比 
 
%%  -------混频器中的低通滤波器设计-------- 
 
  fc=F0;  %截止频率,必须大于w1-w2=Fpian,小于w1+w2=2F0; 
  wc=fc/(Fs/2); % 归一化截止频率 
  M=51;     %滤波器长度,阶数为M-1; 
  W=blackman(M);%Blackman窗 
  hn=fir1(M-1,wc,W);%用Blackman窗生成fir滤波器 
   
%---滤波器的性能指标----- 
  tr_width=11*pi/M;%过度带宽 
  tr_width_fs=(11/M)*(Fs/2);%pi对应Fs/2; 
  delay=(M-1)/2;%群延时 
%   figure;freqz(hn,1,[],Fs); 
 
%%   %%%%%%%%% 锁相环主体程序 %%%%%%%%%%% 
 
%初始化 
Ve_I=zeros(1,num);%I路混频输出 
Ve_Q=zeros(1,num);%Q路混频输出 
fir_out_I=zeros(1,num);%I路FIR滤波器输出 
fir_out_Q=zeros(1,num);%Q路FIR滤波器输出 
loop_out=zeros(1,num);%环路滤波器输出 
phase1=zeros(1,num);%相对于固有振荡频率的输入相位 
phase2 = zeros(1,num);%vco相对于固有振荡频率的输出相位 
 
pd_out=zeros(1,num);%鉴相器输出 
 
vco_out_I=cos(2*pi*Fn.*t+phase2);%Vco输出 
vco_out_Q=sin(2*pi*Fn.*t+phase2);%Vco输出到Q路 
% showfft(vco_out_I,Fs,'vco输出频谱') 
 
Ve_ini_I=vco_out_I.*Signal;%初始混频输出 
% showfft(Ve_ini_I,Fs,'I路初始混频输出信号频谱') 
 
pd_ini_I=filter(hn,1,Ve_ini_I); 
%   showfft(pd_ini_I,Fs,'I路初始鉴相输出信号频谱'); 
 
fd(1)=wn^2*loop_out(1)/(2*pi); 
 
tic; %秒表计时器开始计时 
 
h = waitbar(0,'Please wait...'); 
 
for n=2:num 
    waitbar((n-1)/(num-1),h); 
     
     %混频输出 
     Ve_I(n)=1*vco_out_I(n-1)*Signal(n-1);  %此处可考虑乘以一个增益 
     Ve_Q(n)=1*vco_out_Q(n-1)*Signal(n-1); 
      
     %FIR滤波器输出 
     if n>=M 
            fir_out_I(n)=Ve_I(n:-1:n-M+1)*(hn'); 
            fir_out_Q(n)=Ve_Q(n:-1:n-M+1)*(hn'); 
     else 
            fir_out_I(n)=fliplr(Ve_I(1:1:M))*(hn'); 
            fir_out_Q(n)=fliplr(Ve_Q(1:1:M))*(hn'); 
     end 
      
     %鉴相器输出 
     pd_out(n)=atan(fir_out_I(n)/fir_out_Q(n)); 
 
     %环路滤波器输出 
     loop_out(n)=A*loop_out(n-1)+K1*pd_out(n)-K2*pd_out(n-1); 
      
     %压控振荡器输出 
     phase2(n) = phase2(n-1) + K3*loop_out(n)+K3*loop_out(n-1) ; 
      
     vco_out_I(n)=cos(2*pi*Fn*n*Ts+phase2(n)); 
     vco_out_Q(n)=sin(2*pi*Fn*n*Ts+phase2(n)); 
     
     fd(n)=wn^2*loop_out(n)/(2*pi); 
end 
close(h); 
 
%%   %%%%%%计算多普勒(测速)和测频误差%%%%%%%%%% 
 
%为提高测频精度,对结果进行积分滤波 
jifen=1/3*num;%积分时间为总采样点数的4分之一,即0.25*T秒;合理取积分点数,可以减小跟踪误差 
fd_jifen=wn^2*sum(loop_out(num-jifen:num))/(2*pi*jifen); 
wucha=abs(fd_jifen-Fpian);  
 
%%  %%%%%%%%%%绘图显示结果%%%%%%%%%%%% 
 
% figure;plot(t,Signal); 
% xlabel('t (s)');ylabel('幅度');title('输入信号');grid on; 
 
% figure;plot(t,loop_out); 
% xlabel('t (s)'); ylabel('幅度');title('环路滤波器输出'); 
% legend(['多普勒频偏',num2str(fd_jifen),'Hz, ','跟踪误差',num2str(wucha),'Hz']);grid on; 
 
figure;subplot(211);plot(t,pd_out); 
xlabel('t (s)');ylabel('弧度(rad)'); 
title(char(['鉴相器输出的相位差(S/N=',num2str(snri),'dB 多普勒频偏=',num2str(Fpian),')'],['跟踪误差为',num2str(wucha),'Hz'])); 
legend(['f0=',num2str(F0*1e-3),'KHz fs=',num2str(Fs*1e-3),'KHz BL=',num2str(BL),'Hz']);grid on; 
subplot(212);plot(t,fd); 
xlabel('t (s)');ylabel('频率(Hz)'); 
title(char(['vco偏置频率(S/N=',num2str(snri),'dB 多普勒频偏=',num2str(Fpian),')'],['跟踪误差为',num2str(wucha),'Hz'])); 
legend(['f0=',num2str(F0*1e-3),'KHz fs=',num2str(Fs*1e-3),'KHz BL=',num2str(BL),'Hz']);grid on; 
 
% figure;plot(t,vco_out); 
% xlabel('t (s)');ylabel('幅度');title('vco输出信号');grid on; 
 
% figure;plot(t,phase2); 
% xlabel('t (s)');ylabel('弧度(rad)');title('vco输出相位');grid on; 
 
% %画出输入信号与vco输出信号的频谱图以作比较 
% showfft(vco_out_I,Fs,'vco输出信号频谱'); 
% showfft(Ve_I,Fs,'混频输出信号频谱');