cpafc.m

通信中常用的叉积自动频率频率控制算法

clc;
%参数初始化
N=2048
sum=0;
f0=500;
fs=8000;
phi(1)=0;
temp3=0;
phi(1) = 0;
temp1 = 0;
temp_pre1 = 0;
temp2 = 0;
temp3 = 0;

T=1/fs;
coff=0.0001
BL =coff*fs  %Arbitrarily set
BL=1
BLf=10
Wn=1.89*BL
Wnf=4*BLf
c1=1.414*Wn*T;
c2=Wn*Wn*T
d1=Wnf
zeta = 0.707;  %set in Spec
k1 =1.00;%mean(abs(input.^2))
alpha =2*zeta*2*BL/(zeta+1/4/zeta)/fs/k1 %
beta =((2*BL/(zeta+1/4/zeta)/fs))^2/k1%
%打开文件
% fid=fopen('5751.wav');
% in=fread(fid,44,'short');
% in=fread(fid,N,'short')/32767;
% in=fread(fid,N,'short')/32767;
% in=fread(fid,N,'short')/32767;
% in=fread(fid,N,'short')/32767;
% in=fread(fid,N,'short')/32767;
% in=fread(fid,N,'short')/32767;
% in=fread(fid,N,'short')/32767;
 fid=fopen('am.dat');
 in=fread(fid,N,'float');
 
%低通滤波器设计
a=[1 0];
fcut=[100  300];
dev=[0.1 0.01]
[n,Wn,beta_lp,ftype]=kaiserord(fcut,a,dev,fs);
n=n+rem(n,2);
hd=fir1(n,Wn,ftype,kaiser(n+1,beta_lp),'scale');
%滤波器系数倒置
F_Len=n+1;
for i=1:F_Len
    hds(i)=hd(F_Len-i+1);
end
for i=1:N
%     I(i)=0;
%     Q(i)=0;
    fer(i)=0;
    q(i)=0;
    yyI(i)=0;
    yyQ(i)=0;
end

for i=1:N
    phi(i)=temp3;
    I(i)=in(i)*cos(phi(i));
    Q(i)=in(i)*sin(phi(i));
    %滤波
     if(i>F_Len)
        for j=1:F_Len
         nI(j)=I(i-j+1);
         nQ(j)=Q(i-j+1);
         nin(j)=in(i-j+1);
        end
        yyI(i)=dot(nI,hds);
        yyQ(i)=dot(nQ,hds);
        sig1=-sign(yyI(i-1)*yyI(i)+yyQ(i)*yyQ(i-1));
        fer(i)=(yyI(i-1)*yyQ(i)-yyI(i)*yyQ(i-1))*sig1*1; %频差
     end
%     temp1=temp_pre1+fer(i)*beta;
%     temp2=alpha*fer(i)+temp1;
     f0=f0+d1*fer(i);
     temp3=2*pi*f0*T+phi(i);
%     temp_pre1=temp1;
%     temp2=freq(i)+c2*e(i)+d1*f(i);
%     temp3=2*pi*f0*T+phi(i)+c1*e(i);
%     
    
    
end

f=(0:(N/2-1));
for i=1:N/2-1;
    f(i)=f(i)*fs/N;
end
%输入信号谱
y=fft(in,N);
pyy=y.*conj(y);
subplot(4,1,1);
plot(f,pyy(1:N/2))
title('输入信号频谱')
%I路信号频谱
y=fft(I,N);
pyy=y.*conj(y);
subplot(4,1,2);
plot(f,pyy(1:N/2))
title('I 路信号频谱')
%I路信号频谱
y=fft(Q,N);
pyy=y.*conj(y);
subplot(4,1,3);
plot(f,pyy(1:N/2))
title('Q 路信号频谱')
%滤波器谱
y=fft(hd,N);
pyy=y.*conj(y);
subplot(4,1,4);
plot(f,pyy(1:N/2))
title('滤波器频谱')
figure(2)
%输入   滤波谱
y=fft(yyI,N);
pyy=y.*conj(y);
subplot(3,1,1);
plot(f,pyy(1:N/2))
title('输入滤波谱')
%I路    滤波谱
y=fft(yyI,N);
pyy=y.*conj(y);
subplot(3,1,2);
plot(f,pyy(1:N/2))
title('I 滤波谱')
%Q路    滤波谱
y=fft(yyQ,N);
pyy=y.*conj(y);
subplot(3,1,3);
plot(f,pyy(1:N/2))
title('Q 滤波谱')
%输入信号波形
figure(3)
subplot(3,1,1);
plot(in)
title('原始信号')
xlabel('time')

subplot(3,1,2);
plot(fer)
title('频差')
subplot(3,1,3);
plot(yyI)
title('输出')