Whoru 发表于 2009-6-23 16:33  振动论坛贴子主题 <求教：apFFT频谱校正问题！>、

我想写一个通用的apFFT时移相位差法校正程序，也就是对任意的信号，可以采用apFFT作频谱校正，主要参考论坛里面贴出的程序，可是得到的结果非常不理想，请问高手是什么问题？

http://www.chinavib.com/forum/viewthread.php?tid=83661&extra=page%3D3%26amp%3Bfilter%3Dtype%26amp%3Btypeid%3D186

 

<数字信号全相位谱分析与滤波技术>一书附录二和附录三的2个程序都是用於理解校正原理的程序, 不是实用程序, 所以

1)采样频率选为N, 即fft的频率分辨率为1Hz/格, 不作频率转换使初学者搞不清楚

2)不用判断语句, 所以出现round 和floor语句

3)不加噪音

4)不加自动搜索振幅峰值语句

5)加已知余弦信号, 这样可以知道校正是否正确

 

(1)自动搜索峰值的apfft/apfft校正程序如下(在楼主Whoru 的程序上修改而成)

function [fffa,AA,PH]=apFFT_Corret(k0,F2,F1,N);
clc;clear;clf;close all;
Fs=4096;
N=8192;
t=-N+1:N*2-1;
f=[64.2,126.4,131.7,258.5,382.8];
A=[5 4 3 2 1];
Ph=[40,20,30,70,60];
CNum=length(f);
s=zeros(1,length(t));
for i=1:CNum
myt=A(i)*cos(2*pi*t*f(i)/Fs+Ph(i)*pi/180);
s=s+myt;
end
s=s(1:3*N-1);
win=sin(pi*(1/2:N-1/2)/N).^2;
win2=conv(win,win);
win2=win2/sum(win2);
s1=s(1:2*N-1);
y1=s1.*win2;
y1a=y1(N:end)+[0,y1(1:N-1)];
F1=fft(y1a,N);
s2=s(1+N:3*N-1);
y2=s2.*win2;
y2a=y2(N:end)+[0,y2(1:N-1)];
F2=fft(y2a,N);
a1=abs(F1);
L=fix(N/2)+1;
f=(0:L)*Fs/N;
for i=1:CNum
[A,Ind(i)]=max(a1(1:L));
k=Ind(i);
peak_F2(i)=F2(k);
peak_F1(i)=F1(k);
StartP=Ind(i)-2;
StopP=Ind(i)+2;
if Ind(i)-2<1
StartP=1;
end
if Ind(i)+2>L
StopP=L;
end
Indx=StartP:StopP;
a1(Indx)=zeros(1,length(Indx));
end
num=length(Ind); 
for m=1:num;       
[fff(m) ,AA(m) ,PH(m) ]= apFFT_Corret(Ind(m),peak_F2(m),peak_F1(m),N);
end
disp('校正的频率')
fm=[fff]*Fs/N
disp('校正的相位')
PH=mod(PH*180/pi,360)
disp('校正的振幅')
AA
function [fff ,AA ,PH ]= apFFT_Corret(k0,F2,F1,N)  
dph=angle(F2)-angle(F1);
dph=mod(dph,2*pi);
if dph>pi
dph=dph-2*pi;
elseif dph<-pi
dph=dph+2*pi;
end   
dk=dph/pi/2;
AA=2*((1-dk.*dk)./sinc(dk))^2.*abs(F1);

fff=(k0+dk-1);
PH=angle(F1);

运行结果:

fm = 64.2                           126.399999995859      131.700000003699      258.5                       382.8
PH = 39.9999999999998     20.0000005798297      29.9999994821281     69.9999999999999  60.0000000000007
AA = 5.0000000000001         4.00000024675366      3.00000017716603     2                               1.00000000000003

 

(2)自动搜索峰值的fft/apfft校正程序如下

function [fffa,AA,PH]=apFFT_Corret(k0,F2,F1,N);
clc;clear;clf;close all;
N=8192;Fs=4096;
t=-N+1:N-0;
f=[64.2,126.4,131.7,258.5,382.8];
A=[5 4 3 2 1];
Ph=[40,20,30,70,60];
CNum=length(f);
s=zeros(1,length(t));
for i=1:CNum
myt=A(i)*cos(2*pi*t*f(i)/Fs+Ph(i)*pi/180);
s=s+myt;
end
win=hanning(N)';
win1=win/sum(win);
win2=conv(win,win);
win2=win2/sum(win2);
s1=s(N:2*N-1);
y1=s1.*win1;
F1=fft(y1,N);
s2=s(1:2*N-1);
y2=s2.*win2;
y2a=y2(N:end)+[0,y2(1:N-1)];
F2=fft(y2a,N);
a1=abs(F1);
L=fix(N/2)+1;
f=(0:L)*Fs/N;
for i=1:CNum
[A,Ind(i)]=max(a1(1:L));
k=Ind(i);
peak_F2(i)=F2(k);
peak_F1(i)=F1(k);
StartP=Ind(i)-2;
StopP=Ind(i)+2;
if Ind(i)-2<1
StartP=1;
end
if Ind(i)+2>L
StopP=L;
end
Indx=StartP:StopP;
a1(Indx)=zeros(1,length(Indx));
end
num=length(Ind); 
for m=1:num;       
[fff(m) ,AA(m) ,PH(m) ]= apFFT_Corret(Ind(m),peak_F2(m),peak_F1(m),N);
end
disp('校正的频率')
fm=[fff]*Fs/N
disp('校正的相位')
PH=mod(PH*180/pi,360)
disp('校正的振幅')
AA
function [fff ,AA ,PH ]= apFFT_Corret(k0,F2,F1,N)  
dph=angle(F1)-angle(F2);
dph=mod(dph,2*pi);
if dph>pi
dph=dph-2*pi;
elseif dph<-pi
dph=dph+2*pi;
end   
dk=dph/(pi*(N-1)/N);
AA=abs(F1).^2/abs(F2)*2;
fff=(k0+dk-1);
PH=angle(F2);

运行结果:

fm = 64.1999999907037   126.399970553215   131.700037295151       258.500000006326       382.800000001305
PH =39.9999999999998     20.0000005809033    29.9999994769034     69.9999999999999    60.0000000000007
AA = 4.99999886668805    4.00077883432835    3.00074152949353    2.00000011889194    0.999999966045438

 

(3)一位时移4cfft/4cfft校正算法(下程序中的N=8196/4,以保证输入数据不增多)

4阶循环移位窗fft(4cfft)的相位精度比apfft高,用4cfft构成一位时移相位差法效果好於用apfft构成的一位时移相位差法

function [fffa,AA,PH]=apFFT_Corret(k0,F2,F1,N);
clc;clear;clf;close all;
Fs=4096;
N=8192/4;
t=1:N*4-2;
f=[64.2,126.4,131.7,258.5,382.8];
Ph=[40,20,30,70,60];
A=[11 7 5 3 1]*1;
CNum=length(f);
s=zeros(1,length(t));
for i=1:CNum
myt=A(i)*cos(2*pi*t*f(i)/Fs+Ph(i)*pi/180);
s=s+myt;
end
s=s(1:4*N-2);
win=sin(pi*(1/2:N-1/2)/N).^2;
win2=conv(win,win);
win4=conv(win2,win2);
      s2=s(1:4*N-3);
      s22=s2.*win4;
      sb=[0 0 s22(1:N-2)]+s22(N-1:2*N-2)+s22(2*N-1:3*N-2)+[s22(3*N-1:4*N-3) 0];
      F1=fft(sb,N);
      a1=abs(F1);
s2=s(2:4*N-2);
s22=s2.*win4;
 sb=[0 0 s22(1:N-2)]+s22(N-1:2*N-2)+s22(2*N-1:3*N-2)+[s22(3*N-1:4*N-3) 0];
F2=fft(sb,N);
a1=abs(F1);
L=fix(N/2)+1;
for i=1:CNum
[A1,Ind(i)]=max(a1(1:L));
k=Ind(i);
peak_F2(i)=F2(k);
peak_F1(i)=F1(k);
StartP=Ind(i)-2;
StopP=Ind(i)+2;
if Ind(i)-2<1
StartP=1;
end
if Ind(i)+2>L
StopP=L;
end
Indx=StartP:StopP;
a1(Indx)=zeros(1,length(Indx));
end
num=length(Ind); 
for m=1:num;       
[fff(m) ,AA(m) ,PH(m) ]= apFFT_Corret(Ind(m),peak_F2(m),peak_F1(m),N);
end
fm=[fff]*Fs/N
fm-f
PH=mod(PH*180/pi,360)
PH-Ph
AA
(AA-A)./AA
function [fff ,AA ,PH ]= apFFT_Corret(k0,F2,F1,N)  
dph=mod(phase(F2)-phase(F1),2*pi);
fff=dph*N/2/pi;
PH=angle(F1)-(2*N-1)*dph;
dk=mod(fff,1);
if dk>0.5;
dk=dk-1;
elseif dk<0.5;
dk=dk;
end
AA=2*((1-dk*dk)/sinc(dk))^4*abs(F1)/N/N/N/N*4*4*1;

运行结果:

fm =64.1999999999999   126.3999999949762   131.7000000333688   258.5000000000001   382.8000000000001
PH =40.000000000025466  20.000001911503205  29.99998730354855  69.999999999956344  59.999999999970896
AA =10.999999999999625   6.999999975178254   4.99999990401475  3.000000000000063   0.999999999999961

        apFFT的初始相位是准确的，相位差也是准确的，所以频率可以准确校正,  进而频偏也是准确的.  所以振幅也可以准确校正. 无噪时apFFT的相位,频率,振幅校正精度均达10(-10).