qiugonglvcha.m

解调时考查频率偏移量的一种鉴频算法 采用了双滤波器结构

function [deltP]=qiugonglvcha(of,Nt)
figure(1);
fs=10000;Ts=1/fs;                  %将时间参数设置与主函数一致
t=-10:Ts:10;
k1=(t+10)/10;


[b,a]=fir2(100,[0 0.1 0.1 1],[1 1 0 0]);     %产生滤波器，并滤得基带的窄带信号A(t)
At=filter(b,a,Nt);
Ak=fft(At);
subplot(3,3,2);plot(k1,abs(Ak)); title('A(k)');xlabel('k');

fc=3000;                            %设置调制载波程序
Bt=cos(2*pi*(fc+of)*t);        
Ft=At.*Bt;                          %调制得到调制信号F(t)
Fk=fft(Ft);
subplot(3,3,3);plot(k1,abs(Fk));title('F(k)');xlabel('k');

I=Ft.*cos(2*pi*fc*t);               %对载波信号解调分别得到其实部和虚部成分
Q=Ft*(-1).*sin(2*pi*fc*t);
Ik=fft(I);
subplot(3,3,4);plot(k1,abs(Ik));title('I(k)');xlabel('k');
Qk=fft(Q);
subplot(3,3,5);plot(k1,abs(Qk));title('Q(k)');xlabel('k');
Ct=I+j*Q;                           %解调后得到复信号C(t)
Ck=fft(Ct);
subplot(3,3,6);plot(k1,abs(Ck));title('C(k)');xlabel('k');


h1=fir2(length(t)-1,[0 0.05 0.05 1],[1 1 0 0]);      %产生长度与解调信号相同的窄带滤波器h1(n1)
n1=0:(length(h1)-1);
h1=h1.*exp(j*2*pi*0.1*(fs/2)*Ts*n1);     %令滤波器中心频率由0转到f1=0.1*(fs/2),其对应的移位因子是exp(j*w*n1)=
                                         %exp(j*Ω*Ts*n1)=exp(j*(2*pi*f1)*Ts*n1)=exp(j*(2*pi*0.1*(fs/2))*Ts*n1)
H1=fft(h1);                                   %复数滤波器冲击响应的FFT
k11=n1/(length(h1)/2);                        %指定该FFT在频域的范围k11，并将其归一化到区间[0，2]
subplot(3,3,7);plot(k11,abs(H1));title('H1(k)');xlabel('k');

S1k=Ck.*H1;                                   %在频率中进行滤波并得到上支路信号总能量         
P1xulie=abs(S1k).^2;
P1=sum(P1xulie);


n2=n1;                                        %设置中心频率向负方向偏移的滤波器h2(n2),且与h1(n1)的中心频率对称
h2=h1.*exp(j*2*pi*(-0.1)*n2);
k22=n1/(length(h2)/2); 
H2=fft(h2);
subplot(3,3,8);plot(k22,abs(H2));title('H2(k)');xlabel('k');

S2k=Ck.*H2;
P2xulie=abs(S2k).^2;
P2=sum(P2xulie);

if (of<=100)&(of>=-100);     %产生频偏在－100，0，100处两复数滤波器的输出信号的频谱
figure(2)
subplot(2,3,(of+200)/100);plot(k11,abs(S1k));
subplot(2,3,(of+200)/100+3);plot(k22,abs(S2k));
end


deltP=P1-P2                                 %求能量差deltp