zhwang554的个人博客分享 http://blog.sciencenet.cn/u/zhwang554

博文

4阶移位循环卷积窗fft

已有 8935 次阅读 2011-7-14 19:56 |个人分类: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]      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)为取主值区间

 

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

    所以apfft是加2阶移位循环卷积窗fft  

 

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

 

   三角窗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次方

 

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

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

   所以移位循环卷积窗fft在输入数据和性能相同时, 其计算量小於线性卷积窗fft, 相位谱不须校正 

   

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

 

                               图2    3hanning移位循环卷积窗fft

                           3    3hanning线性卷积窗fft

      23hanning移位循环卷积窗fft, 输入数据3N-2, ], N=3移位相加后, 作NFFT

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

 

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

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

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

    高阶线性卷积窗fft和高阶移位循环卷积窗fft比较, 两者输入输入数据相同, 振幅谱连续曲线完全相同, 计算精度相同线性卷积窗fft作p*NFFT, 循环卷积窗fftNFFT,循环卷积窗fft计算量小於线性卷积窗fft

      高构p阶移位循环卷积窗只要有足够多的数据,其泄漏p次方下降, 而FFT阶数并不随数据加长而增加, 是功率谱估计的工具

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

 

                                             图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);  

       4阶移位循环卷积窗fft的相位谱能保持相位不须校正和水平相位特性.这点一开始不信,因为4阶循环卷积窗fft是4段N长的数据相加而成, 最后形成的N个数据中的每一个数据包含3-4个原始采样数据, 相位谱测的是多个原始采样相位值的加权和.
      4阶移位循环卷积窗fft中相位谱测的是4N-3长数组中中间采样点(笫2N-1点)的相位
     在4阶循环卷积窗fft中数据形成十分重要,以下二条差一点都不行;
1 t=-2*N+2:2*N-2,确保中间值为初相位
2 4阶循环卷积窗N为周期移位相加中有4种形成方式:
  s1=[s(1:N)+s(N+1:2*N)+s(2*N+1:3*N)+[s(3*N+1:4*N-3) 0 0 0];
  s2=[0 s(1:N-1)]+s(N:2*N-1)+s(2*N:3*N-1)+[s(3*N:4*N-3) 0 0];
  s3=[0 0 s(1:N-2)]+s(N-1:2*N-2)+s(2*N-1:3*N-2)+[s(3*N-1:4*N-3) 0];
  s4=[0 0 0 s(1:N-3)]+s(N-2:2*N-3)+s(2*N-2:3*N-3)+[s(3*N-2:4*N-3)];
     只有s3保持相位不须校正和水平相位特性,即笫一段前补2个零,笫四段后补1个零另外3种都是左右倾斜线,数值也不对
     s3的笫一个象素是s(N-1)+s(2N-1)+s(3N-1)三个采样的相位加权和, 但s(N-1和s(3N-1)的相位相对中问点s(2N-1)正好相反, s(2N-1)和s(3N-1)的窗加权值又正好相同,所以笫一个象素相位值只由s(N-1)采样的相位决定,从而保持初相位不须校正和水平相位特性

      2,4,6阶等偶数阶移位循环卷积窗fft的相位谱均和apfft一样, 无噪时都具有水平相位特性,且保持相位不须校正.由於水平相位特性是峰值泄漏产生的, 有噪时,2cfft的峰值泄漏仍较大, 仍具有水平相位特性, 而4,6,8阶峰值泄漏太小,水平相位特性消失. 由於水平相位特性可用於判断频谱峰值是否只是一个频率成分,这个有用性质只有apfft才具有.

      1,3,5阶等奇数阶移位循环卷积窗fft的相位谱均和传统1阶fft一样,相位须校正

      4阶移位循环卷积窗fft(4cfft)的相位校正精度高, 可用於一位时移的4cfft/4cfft校正法, 虽相位精度比N位时移4cfft/4cfft校正法要低, 但总的相位精度也较高, 但输入数据少了, 所以一位时移的4cfft/4cfft校正法是一个较好的选择.  如N=256的一位时移4cfft法只须1021+1=1023个输入. 时移法须二次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阶或6阶)移位循环卷积窗fft.能确定相隔一个频率点的二个相邻频率成分的频率相位和振幅. 如

图7所示的相隔一个频率点的6个相邻频率成分, 这是不能再靠近的密集频谱了, 实现了密集频谱分析和校正.

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

 

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

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

         apfft/apfft_1   N=32  fs=32         cfft/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窗的log10频谱图

 

 

求相位差用APFFT,对频谱间隔2个以上频率分辨率信号校正精度高.

                  4cFFT,对频谱间隔1个以上频率分辨率信号校正精度高.

   

 

 

 

 

 

 



https://blog.sciencenet.cn/blog-212920-464833.html

上一篇:全相位非整周期采样三角函数的正交性和整周期性--apDIQ
下一篇:4cfft插值公式证明1
收藏 IP: 74.105.210.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-3-28 21:43

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部