日志
4阶移位循环卷积窗fft
||
卷积有两种: 线性卷积和循环卷积。卷积窗也有两种: 线性卷积窗和循环卷积窗。卷积定理也有两种: 线性卷积定理和循环卷积定理。加卷积窗可使插值公式平方，加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]
F1长3, 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]
F1长3, F4长3, 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)为取主值区间
还可构成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] 。
F1的3方为F5 . 即3阶循环卷积窗的插值公式是加三角窗fft插值公式的3次方。
但是这里不是把循环卷积窗直接加到信号上,而是由线性卷积移位形成的循环卷积窗.循环卷积窗直接加到0:N-1的信号上,在整数倍取样时也有平方振幅谱和水平相位谱,但在非整数倍取样时就没有平方振幅谱和水平相位谱。而移位循环卷积窗fft,在整数倍取样时有平方振幅谱和水平相位谱,在非整数倍取样时也有平方振幅谱和水平相位谱。
所以apfft是加2阶移位循环卷积窗fft。是移位循环卷积窗的最低阶fft。 可以构成3阶,4阶等高阶移位循环卷积窗fft。
图1是N=64的 1 2 4阶hanning移位循环卷积窗fft的对数振幅谱和相位谱,.其输入数据分别是64 127, 253个,均作64阶FFT。.偶数阶移位循环卷积窗fft相位谱呈水平特性, 且不用校正。
相比N=64的 1 2 4 阶hanning线性卷积窗fft, 其输入数据分别是64 127, 253个,其对数振幅谱类同图一,但分别是64,128,256阶FFT,
图1 1 2 4 阶hanning移位循环卷积窗的对数振幅谱
图2 3阶hanning移位循环卷积窗fft
图3 3阶hanning线性卷积窗fft
图2是3阶hanning移位循环卷积窗fft, 输入数据3N-2, 经N为周期移位相加后, 作N阶FFT。
图3是3阶hanning线性卷积窗fft. 输入数据3N-2,作3N阶FFT。
图2和图3输入数据相同, 性能类同, 移位循环卷积窗fft计算量小於线性卷积窗fft, (3阶hanning线性卷积窗没有快速运算, 4阶hanning线性卷积窗有快速运算, 4阶图类同)。
图4 3阶移位循环卷积窗fft和3阶线性卷积窗fft的对数振幅谱
图4 是3阶移位循环卷积窗fft(图3b)和3阶线性卷积窗fft(图3a)的对数振幅谱, 两者的连续曲线完全相同, 只是3阶线性卷积窗fft的点数比3阶循环卷积窗fft的点数密3倍。所以3阶线性卷积窗fft虽用3N阶FFT得3N阶的频谱, 但其实际频谱分辨力并没有增加, 这和fft中补零的情况一样, 只是频谱变密。
图5是4阶hanning移位循环卷积窗fft方框图,图6是N=256.fs=3000时4阶hanning移位循环卷积窗fft振幅谱和相位谱,其比值法程序见下附.
4阶hanning移位循环卷积窗fft的相位谱仍呈水平特性,其值不须校正.
2,4,6,8阶等偶数阶移位循环卷积窗fft的相位谱均和apfft一样, 具有水平相位特性,且保持相位不须校正
图5 4阶hanning移位循环卷积窗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
