我要加入 登录
声振论坛 返回首页

zhwang554的个人空间 http://home.vibunion.com/?62061 [收藏] [复制] [分享] [RSS]

日志

4阶移位循环卷积窗fft

已有 1898 次阅读2012-6-3 00:32 |个人分类:apfft|

       

 

卷积有两种: 线性卷积和循环卷积卷积窗也有两种: 线性卷积窗和循环卷积窗。卷积定理也有两种: 线性卷积定理和循环卷积定理卷积窗可使插值公式平方,加p阶卷积窗可使插值公式p次方, 泄漏成倍减小

       DFT是离散Fourier変换, 适用於循环卷积定理, 而不是线性卷积定理, 所以DFT中直接加线性卷积窗引起主辨宽度变宽, 降低了分辨力。应该在DFT中加循环卷积窗。全相位中的移位相加,将长2N-1线性卷积窗变成长N循环卷积窗,实现了在DFT中加卷积窗, 使插值公式的平方,而窗主辨宽度不変。

       直观来看, N阶窗的频谱长N, 两个N阶的线性卷积窗的频谱长2N-1, 所以2N-1阶的线性卷积窗的频谱不可能是N阶窗的频谱的平方, 而是补零后的窗的频谱的平方。

         两个N阶的循环卷积窗的频谱长仍是N, 所以N阶的循环卷积窗的频谱才是N阶窗的频谱的平方。

    

     用一个简单例子说明线性卷积窗和循环卷积窗的区別和关联:

N=3三角窗                                     w1=[1 2 1]               DFT振幅为      F1=[4  1  1]

N=3三角窗线性卷积窗    w2=[1 4 6 4 1 ]       DFT振幅为      F2=[16     6.8541   0.1459     0.1459     6.8541]

          F13,  F2长5    F2不可能是F1的平方

补零三角窗                                         w3=[1 2 1 0 0 0]      DFT振幅为     F3= [4   2.616   0.3819   0.3819    2.616]

         F3长5 F2长5,    F3的平方为F2     线性卷积窗频谱是补零三角窗频谱的平方

三角窗循环卷积窗为                         w4=[6 5 5]                DFT振幅为      F4=[16  1  1]

        F13,  F43,   F1的平万为F4      循环卷积窗频谱是三角窗频谱的平方

       线性卷积窗和循环卷积窗的关联如下:

    线性卷积窗为w2=[1 4 6 4 1], 经周期移拓后 (周期移拓是把原序列平移某个周期T的整数倍后再全部加起来所产生的新序列),即笫1位和第4位相加, 2位和第5位相加), 即得循环卷积窗[6 5 5]

线性卷积和循环卷积之间的关系数学上描述如下

x1(n)x2(n)的线性卷积为yl(n), x1(n)x2(n)的循环卷积yc(n)yl(n)的周期延拓,下式中GN(n)为取主值区间

上式中GN(n)为取主值区间

         还可构成3,4阶等高阶循环卷积窗fft

         N=3三角窗3阶线性卷积窗为w5=[1 2 1]*[1 2 1]*[ 1 2 1]=[1 6 15 20 15 6 1]

          3线性卷积窗w5=[1 6 15 20 15 6 1 ], N=3周期移拓(即笫1,4位和第7位相加, 2位和第5位相加, 3位和第8位相加) 得三角窗3 阶循环卷积窗w6=[22 21 21] DFT振幅为 F5-=[64 1 1]

           F13方为F5 . 即3阶循环卷积窗的插值公式是加三角窗fft插值公式的3次方

            但是这里不是把循环卷积窗直接加到信号上,而是由线性卷积移位形成的循环卷积窗.循环卷积窗直接加到0:N-1的信号上,在整数倍取样时也有平方振幅谱和水平相位谱,但在非整数倍取样时就没有平方振幅谱和水平相位谱而移位循环卷积窗fft,在整数倍取样时有平方振幅谱和水平相位谱,在非整数倍取样时也有平方振幅谱和水平相位谱

          所以apfft是加2阶移位循环卷积窗fft。是移位循环卷积窗的最低阶fft 可以构成3,4阶等高阶移位循环卷积窗fft

         图1是N=64 1 2 4hanning移位循环卷积窗fft的对数振幅谱和相位,.其输入数据分别是64 127, 253,均作64FFT.偶数阶移位循环卷积窗fft相位谱呈水平特性, 且不用校正

相比N=64 1 2 4 hanning线性卷积窗fft, 其输入数据分别是64 127, 253,其对数振幅谱类同图一,但分别是64,128,256FFT,

                           图1         1 2 4 hanning移位循环卷积窗的对数振幅谱

                                                                 2    3hanning移位循环卷积窗fft

                                                  3             3hanning线性卷积窗fft

     图23hanning移位循环卷积窗fft, 输入数据3N-2, N为周期移位相加, NFFT

     图33hanning线性卷积窗fft. 输入数据3N-2,3NFFT

     图2和图3输入数据相同, 性能类同, 移位循环卷积窗fft计算量小於线性卷积窗fft, (3hanning线性卷积窗没有快速运算, 4hanning线性卷积窗有快速运算, 4阶图类同)

                      4      3阶移位循环卷积窗fft3阶线性卷积窗fft的对数振幅谱

      图4 3阶移位循环卷积窗fft(3b)3阶线性卷积窗fft(3a)的对数振幅谱, 两者的连续曲线完全相同, 只是3阶线性卷积窗fft的点数比3阶循环卷积窗fft的点数密3。所以3线性卷积窗fft虽用3NFFT3N阶的频谱, 但其实际频谱分辨力并没有增加, 这和fft中补零的情况一样, 只是频谱变密

     5是4hanning移位循环卷积窗fft方框图,图6是N=256.fs=3000时4hanning移位循环卷积窗fft振幅谱和相位谱,其比值法程序见下附.

     4阶hanning移位循环卷积窗fft的相位谱仍呈水平特性,其值不须校正.

      2,4,6,8阶等偶数阶移位循环卷积窗fft的相位谱均和apfft一样, 具有水平相位特性,且保持相位不须校正

     

 

 

                         图5                4hanning移位循环卷积窗fft

                                图6     4cfft对数振幅谱相位谱  

4阶移位循环卷积窗fft比值法程序

function XfCorrect=SpectrumCorrect(N,xf,CorrectNum,fs)
close all;clear all;clc 
N=256;fs=3000;
t=(-2*N+2:N*2-2)/fs;f=50.0;
s=240*cos(2*pi*f*t+10*pi/180)+0.2*cos(2*pi*2*f*t+20*pi/180)+9*cos(2*pi*3*f*t+30*pi/180)+0.1*cos(2*pi*4*f*t+40*pi/180)+39*cos(2*pi*5*f*t+50*pi/180)+0.02*cos(2*pi*6*f*t+60*pi/180)+21*cos(2*pi*7*f*t+70*pi/180)+0.05*cos(2*pi*8*f*t+80*pi/180)+0.4*cos(2*pi*9*f*t+90*pi/180)+0.08*cos(2*pi*10*f*t+100*pi/180)+0.6*cos(2*pi*11*f*t+110*pi/180);
xx=[f 240 10 ;5*f 39 50; 7*f 21 70;3*f 9  30; 11*f 0.6 110  ;9*f 0.4 90;2*f 0.2 20; 4*f 0.1  40; 10*f 0.08 100;8*f 0.05 80  ;6*f 0.02 60];
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];
      xf=fft(sb,N);
      a1=abs(xf);p1=mod(angle(xf)*180/pi,360);
      p1(5)
      p1(14)
      xf=xf(1:N/2);
XfCorrectW_apfft_2N=SpectrumCorrect1(N,xf,11,fs);
XfCorrectW_apfft_2N(:,1)=XfCorrectW_apfft_2N(:,1);
XfCorrectW_apfft_2N
XfCorrectW_apfft_2N-xx

subplot(211),plot(10*log10(a1),'k.-');
grid
title('N=256   the 4 order cycle convolution hanning window amplitude spectrum');
%ylim([-170,0]);
xlim([0,N/2]);
legend('4cfft');
xlabel('f');ylabel('db');
subplot(212),plot(p1,'k.-');
grid
title('N=256 the 4 order cycle convolution hanning window angle spectrum');
ylim([0,150]);
xlim([0,N/2]);
legend('4cfft');
xlabel('f');ylabel('degree');

 function XfCorrect=SpectrumCorrect1(N,xf,CorrectNum,fs)
XfCorrect=zeros(CorrectNum,3);
for i=1:CorrectNum
    A=abs(xf);
    [Amax,b]=max(A);b;
    phmax=angle(xf(b));  
    if A(b-1)>A(b+1);
        a=A(b)/A(b-1);
        df=(2-a^(1/4))/(1+a^(1/4))*1;
        XfCorrect(i,1)=(b-1-df)*fs/N;
        XfCorrect(i,2)=(1-df.^2).^4/(sinc(df).^4)/N/N/N/N*Amax*4*4*2;
        XfCorrect(i,3)=mod(phmax*180/pi,360);
        else A(b-1)<A(b+1);
        a=A(b)/A(b+1);
      df=(2-a^(1/4))/(1+a^(1/4))*1;
        XfCorrect(i,1)=(b-1+df)*fs/N;
        XfCorrect(i,2)=(1-df.^2).^4/(sinc(df).^4)/N/N/N/N*Amax*4*4*2;
        XfCorrect(i,3)=mod(phmax*180/pi,360);
    end
        xf(b-2:b+2)=zeros(1,5);
    end
    XfCorrect(i,3)=mod(XfCorrect(i,3),360);
    XfCorrect(i,3)=XfCorrect(i,3)-(XfCorrect(i,3)>180)*360;

 

      4阶移位循环卷积窗fft(4cfft)的相位校正精度高, 可用於一位时移的4cfft/4cfft校正法, 虽相位精度比N位时移4cfft/4cfft要低, 但总的相位精度也较高, 可N位时移apfft/a[fft的最大缺点输入数据少了, 所以一位时移的4cfft/4cfft校正法是一个较好的选择.  如N=256的4cfft只须1021+1=1022个输入,虽时移位须二次FFT, 但二次实FFT可合成一次复数FFT. 一位时移的4cfft/4cfft校正法程序参见振动论坛日志

http://home.chinavib.com/home.php?mod=space&uid=62061&do=blog&quickforward=1&id=18060

<自动搜索峰值的 apfft/apfft ,fft/apfft 和一位时移4cfft/4cfft校正程序>

            一位时移高阶(4阶或8阶)移位循环卷积窗fft.能确定相隔一个频率点的二个相邻频率成分的频率相位和振幅. 如

图7所示的相隔一个频率点的6个相邻频率成分,

                        图7                  相隔一个频率点的6个相邻频率成分的对数振幅谱

           4cfft的可校正密集频谱,  在极低频率测量中将发挥作用, 它可抑止镜像频率的干扰, 在频谱上只要>0.5的谱峰就可测量, >1的谱峰测量有较大精度(频率振幅相位精度均在10(-5)以止). 下表一显示一位时移apfft/apfft和一位时移4cfft/4cfft在极低频时的误差. 在f=0.775时,4cfft仍可测定.

表一     一位时移apfft/apfft和一位时移4cfft/4cfft在极低频时的误差.

Apfft/apfft_1    N=32   fs=32                             4cfft/4cfft_1  N=16  fs=16

f=1.15   a=1   p=200

-6.2764e-004  2.4195e-001 -2.5040e-005         -2.4300e-007 -7.7446e-006  1.6895e-004

f=0,875  a=1   p=100

5.0070e-004 - 2.1690e-001  -3.9412e-004        -7.9718e-007 -4.4183e-006  5.6486e-004

f=0.775  a=1   p=100

-3.3157e-003  1.0835e+000  4.4078e-003        1.1730e-005 -3.9668e-005 -7.9784e-003

f=0.675  a=1   p=100

-6.7500e-001 -1.0000e+002  4.1600e-002        -2.6027e-005  5.1837e-005  1.3845e-002

                             图8 是1,2,4阶循环卷积矩形窗的频谱图, 由图可见,高阶循环卷积窗频谱的主辨宽度变窄, 在0.707处,1cfft宽约2.7652, 4cfft宽约1.4

 

                                  图8                    1,2,4阶移位循环卷积矩形窗的频谱图

            图9是1,2,4阶移位循环卷积hanning窗的对数频谱图,可见阶数越高, 主辨宽度变窄, 旁辨峰值越低,旁辨渐进衰减率越高

                                  图9                    1,2,4阶移位循环卷积hanning窗的对数频谱图

求相位差那用APFFT没错的!对频谱间隔2个以上频率分辨率信号校正精度非常高!

          4cFFT没错的!对频谱间隔1个以上频率分辨率信号校正精度非常高!

        但是这里不是把循环卷积窗直接加到信号上,而是由线性卷积移位形成的循环卷积窗.无移位循环卷积窗直接加到0:N-1的信号上,在整数倍取样时也有平方振幅谱和水平相位谱(图10(a)(b)),但在非整数倍取样时就没有平方振幅谱和水平相位谱(图11(a)(b))而移位循环卷积窗fft,在整数倍取样时有平方振幅谱和水平相位谱(图11(c)(d)),,在非整数倍取样时也有平方振幅谱和水平相位谱(图11(c)(d))

                                          图10   f=29.0  hanning 移位循环卷积fft和 hanning 非移位循环卷积fft的振幅谱和相位谱

                      图11   f=29.3  hanning 移位循环卷积fft和 hanning 非移位循环卷积fft的振幅谱和相位谱

  全相位fft是加循环卷积窗fft

 

评论 (0 个评论)

facelist doodle 涂鸦板

您需要登录后才可以评论 登录 | 我要加入

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-12-25 15:23 , Processed in 0.045489 second(s), 16 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

返回顶部