|
最近看到网易上的两篇 一箫一剑走天下 关于apfft日志,转贴这儿,并作些讨论
全相位FFT算法的实现 - 一箫一剑走天下的日志 - 网易博客
一箫一剑走天下 |
http://blog.163.com/fdy_001/blog/static/12010203320091022102250697/
全相位FFT算法的实现
2009-11-22 22:22
关键字:apFFT 频谱校正 全相位
传统FFT利用三角函数的正交性,将信号分离出来,从而将时域的信号变换到频域。但是,它有一个很重要的前提:输入的序列必须是周期内等间隔采样的值,这样,FFT计算的结果才是我们想要的。
实际的情况是,很难做到等间隔采样。比如,交流电的频率是变化的,并不是固定的50Hz。如果采用按照50Hz的信号来采样,则计算结果将无法反映原始信号。
为什么会出现上述的情况呢?因为FFT是将截断的序列做周期延拓而得到无限长的序列的。当不是等间隔采样时,做周期延拓后,在首尾相接的地方,就会出现信号的跳变,与原始信号不一致,自然也就得不到想要的结果,在频谱上的表现就是频谱出现泄露。
传统FFT的结果可以通过一些算法实现频谱校正,如全能重心法、比值法等等。但都是基于FFT的结果,精度有限。
全相位FFT算法是天津大学的王兆华和候正信教授提出的,具有初始相位不变和有效防止频谱泄露的特性,在他的相关书籍上做过严格的证明,并有相应的matlab程序。
我的项目中,需要分析50±0.5Hz工频信号,所以,研究apFFT算法,希望解决非同步采样的频谱校正问题。
apFFT的算法是先对数据进行全相位预处理,然后进行传统的FFT运算。FFT运算大家都熟悉,在这里,我说说全相预处理的过程。
全相位预处理过程如下:
1,构成一个N点的汉宁窗。
2,汉宁窗对自己求卷积,得到2N-1点的卷积窗。
3,求2N-1点的卷积窗的和。
4,将卷积窗的每一项除以卷积窗的和,得到2N-1点的归一化卷积窗。
5,将数据的1:2N-1项和归一化卷积窗相乘,得到加窗的2N-1项。
6,将第1项和N+1项,第2项和N+2项 ... 第N-1项和第2N-1项相加,得到经过全相预处理的N点序列。
接下来就只需要进行FFT,就可以得到apFFT的结果。
如果我们观察结果全相位预处理的信号,由原来的带有跳变的信号变为平滑连续的信号,使得周期延拓后不会出现信号的跳变。
close all;clc;clear all;
N=256;
t=-N+1:2*N-1;
f1=39.8;
f2=51.345;
A1=5;
ph1=30;
s=A1*cos(1*(2*pi*t*f1/N/50-ph1*pi/180));%+A1*cos(1*(2*pi*t*f2/N+ph1*pi/180));
win=hanning(N)';
win1=hanning(N)';
win2=conv(win,win1);
win2=win2/sum(win2);
w=pi*2;
s1=s(1:2*N-1);
y1=s1.*win2;
y1a=y1(N:end) + [0 y1(1:N-1)];
Out1=fft(y1a,N);
a1=abs(Out1);
p1=mod(phase(Out1),2*pi);
s2=s(1+N:3*N-1);
y2=s2.*win2;
y2a=y2(N:end) + [0 y2(1:N-1)];
Out2=fft(y2a,N);
a2=abs(Out2);
p2=mod(phase(Out2),2*pi);
g=mod((p2-p1)/pi/2,1);
h=2*pi*g.*(1-g.*g)./sin(pi*g);
aa1=abs((h.^2).*a2)/2;
rr=round(f1);
disp('频率校正值')
fff=floor(f1)+g(rr+1)
disp('振幅校正值')
aaa=aa1(floor(f1)+1);
disp('初相位校正值')
ppp=p1(rr+1)*180/pi;
通过在VC下编码验证,apFFT具有很高精度的初始相位。采用apFFT/apFFT时移相位差法,可以高精度的修正频率偏差和幅值偏差。这个算法对于非同步采样处理的意义,不言而喻。与传统fft频谱校正的主要区别在于,apFFT的初始相位是准确的,相位差也就是准确的,频率也就可以校正准确,所以具有更高的精度。王教授对各种校正方法做过详细的测试和对比,若希望了解相关内容,可查阅《数字信号全相位谱分析与滤波技术》一书。
一箫一剑走天下
全相位FFT在VC下的测试结果和分析
http://blog.163.com/fdy_001/blog/static/12010203320091027111614878/
2009-11-27 23:16
关键字:全相位 FFT 校正精度
通过对apFFT算法在VC下的编码实现,证实了一些问题,也发现了一些问题。现在将我的测试过程和结果做一个叙述,以便为有需要的人提供一点参考。
首先,可以肯定的是,这个算法的意义,非同重大,毋庸置疑。在此,对王兆华教授及其团队致以崇高的敬意。但是,书中的一些数据和说法略失偏颇。我对信号处理知之甚少,工作的需要,做了一些了解,在此发表一点自己的观点和看法,欢迎指正。
我测试的结果有两大发现:1,测试结果显示,apFFT确实有初始相位不变性,精度非常高,具体多高,后面有相关的测试数据。
2,校正的精度并没有作者所说的那么高。具体的校正精度到底怎么样,没有详细的表述。
作者在书中说,apFFT对相位、频率和幅值都可以高精度的校准。从他的测试方式和相关的数据来看,这种说法是不严密的。
理由如下:
1,他的测试都是基于FFT的分辨率为1Hz的条件来测试的,这一点可以从数据抽取函数看出。
2,他的测试只是测试了很多离散的频点,不具有一般性。
首先说一下,我的算法的结果和matlab的结果一致,也就是说我的算法本身是没有问题的。
我最开始试图采用apFFT/apFFT来做校正,测试发现,相位和频率的精度非常高,但是,幅度的校正精度最差达0.38%。由于我的测试系统的精度要求是0.2%,这达不到我的要求。为此,我在振动论坛上发帖,希望得到答案。但是,最近我老是上不了论坛,不知道怎么回事?
于是,我寻求别的解决办法,通过对比FFT/apFFT和apFFT/apFFT的校正公式,发现它们不同。于是,我用FFT/apFFT算法,再做了一次测试,发现幅度的校正精度比apFFT/apFFT的高,满足我的工程需要。至此,算是真正解决了非同步采样的FFT校正算法的问题。
我测试的前提是:测试的信号时直流加单频信号,采样点数均相同,所有的数据类型都是float型。对于多频率的信号,只要分辨率足够高,是完全可以分离的,这一点不用怀疑。我的测试方法和王教授的有些不同,具体如下:
1,使用不同的FFT分辨率做测试。
2,选取多个频点做测试。我的频点选取更具一般性,选取50.0到60.0,频率步进为0.1Hz。由于校正的原理是根据相位差来校正的,而为整数频点时相位差为0,随着频率的增大,相位差增大。也就是说,我选取的频点覆盖率很多的相位差,因此,更具一般性。
测试结果如下:
序号 算法 频率范围 分辨率 采样点数 频率精度 幅值精度 相位精度
1 FFT/APFFT 50:0.1:60 16 256 2.77e-4 2.05e-3 1.67e-6
2 FFT/APFFT 50:0.1:60 10 256 7.27e-5 4.22e-4 2.2e-7
3 FFT/APFFT 50:0.1:60 8 256 2.55e-4 2.05e-4 6.7e-8
4 FFT/APFFT 50:0.1:60 4 256 1.56e-4 2.96e-5 3.3e-7
5 FFT/APFFT 50:0.1:60 1 256 3.87e-5 2.85e-6 2.77e-7
对一组频率做校正后,我取最大偏差来计算校正精度。从数据上看,精度似乎不是随分辨率的降低而降低的。这是因为我取的频点少,数据不能较准确的反映真实的最大偏差。但是,上述结果还是能说明一些问题。
分析测试结果不难发现:
1,随着分辨率的升高,相位、频率和幅值的校正精度也逐渐升高。
2,相同的分辨率,相位的精度最高,频率的精度和幅值的精度差不多。
3,随着分辨率的降低,幅值的精度略低于频率,相差越来越大。
王教授得出精度可以达到10的-12次方的结果的原因是,他选取的频率多接近整数频率,而一般接近整数频率的相差小,所以,误差也就小,精度自然也就高。当然,很多实际的情况是频率偏差较小,多在整数附近,因此,我们有时可以这样认为。
由于时间紧迫,没有做更多频点的分析,同时,我认为我的这个测试结果已经可以说明问题,不必再做更细的分析。
FFT的分辨率为:采样率除以一次FFT计算的点数。比如,采样率是512sps,数据进行256点FFT,那么FFT的分辨率就为2Hz,也就是说FFT变换的结果之间的间隔是2Hz。同时,我们也看到,一秒钟只能采样2组数据,能分析的最大频率为256Hz。因此,要想分析更高的频率,需要更高的采样率,而要想更高的分辨率,必须采样足够长的时间,他们是相互矛盾的。我们在实际工程中,希望尽快的计算得到结果,也希望有足够高的分辨率,只能权衡而确定一个合理的采样率和FFT计算点数。我们的系统的精度要求是0.2%,因此,选择分辨率为8Hz,较为合适。如果你的系统中需要分析两个频点很近的信号,建议保证两根谱线间的间距大于3.
感谢 一箫一剑走天下 对我们工作的支持.在网易评语中写不上,在这儿回答一下.
"全相位FFT算法的实现"一文程序中振幅校正语句有改进之处
这个程序也出於书中附录三,当时为了避免出现判断语句,在频偏>0.5时,不在振幅峰值处(round(f1)+1)取值,而取最大点前面的一个(floor(f1)+1), 有时影响精度. 注意round语句对(f1) 4舍5入,floor语句对(f1) 取整数.
可参见下程序(加空格后), 它始终在峰值处取值, 在频偏>0.5时可提高精度
close all;clc;clear all;
N=256;
t=-N+1:N*2-1; f1=39.8;A1=5;ph1=30;
s=A1*cos(1*(2*pi*t*f1/N+ph1*pi/180));
win=hanning(N)';win1=hann(N)';
win2=conv(win,win1);
win2=win2/sum(win2);
w=pi*2;
s1=s(1:2*N-1);%第1组(2N-1)个数据
y1=s1.*win2;
y1a=y1(N:end)+[0 y1(1:N-1)];
Out1=fft(y1a,N);
a1=abs(Out1);
p1=mod(phase(Out1),2*pi);
s2=s(1+N:3*N-1); %第2组(2N-1)个数据
y2=s2.*win2;
y2a=y2(N:end)+[0 y2(1:N-1)];
Out2=fft(y2a,N);
a2=abs(Out2);
p2=mod(phase(Out2),2*pi);
g=mod((p2-p1)/pi/2,1);
h=2*pi*g.*(1-g.*g)./sin(pi*g);
aa1=abs((h.^2).*a2)/2;
rr=round(f1)
disp('频率校正值');
fff=floor(f1)+g(rr+1);
disp('振幅校正值');
aaa=aa1(floor(f1)+1)
disp('初相位校正值');
ppp=p1(rr+1)*180/pi;
dph=mod(p2-p1,2*pi);
dph=dph(rr+1);
if dph>pi
dph=dph-2*pi;
elseif dph<-pi
dph=dph+2*pi;
end
dk=dph/pi/2;
g=dk;
h=2*pi*g.*(1-g.*g)./sin(pi*g);
AA=abs((h.^2).*a1)/2;
aaa=AA(rr+1)
另外<数字信号全相位识谱分析和滤波技术>一书中p.68指出:"经进一步实验发现,若采用经hanning窗和hann窗卷积而成的窗进行双窗全相位FFT谱分析和校正,在无噪下的幅值校正精度还可提高3个数量级."
利用上程序可测 峰前取值, 峰处取值,改为hann,和hanning卷积窗的作用
测频偏0.1-0.9时apfft/apfft校正精度, N=256, 取样频率=N, f1=50.1-50.9,如表一
第二列为加双hanning卷积窗, 峰前取值, 频偏0.4-0.9误差大
第三列为加双hanning卷积窗, 峰处取值, 频偏0.6-0.9有改善
第四列为加hann hanning卷积窗, 峰前取值, 精度好,比第二列提高3个数量级
第五列为加hann hanning卷积窗, 峰处取值, 精度最好,比第三列提高3个数量级
所以apfft/apfft振幅校正采用峰处取值, 频偏0.4-0.9精度有改善,加hann hanning卷积窗, 精度可提高3个数量级,
所以apfft/apfft振幅校正采用hann hanning卷积窗,峰处取值
理论上应加双hanning卷积窗, 为什么加hann hanning卷积窗精度好,不清楚,是偶然发现的
表一 apfft/apfft振幅校正精度
|
峰前取值, 双hanning卷积窗 |
峰处取值, 双hanning卷积窗 |
峰前取值, hann,hannjng卷积窗 |
峰处取值, hann,hannjng卷积窗 |
0.1 |
-1.0109e-004 |
-1.0109e-004 |
-1.9754e-007 |
-1.9754e-007 |
0.2 |
-4.0588e-004 |
-4.0588e-004 |
-7.9939e-007 |
7.9939e-007 |
0.3 |
-9.1897e-004 |
-9.1897e-004 |
-1.8342e-006 |
-1.8342e-006 |
0.4 |
-0.0016 |
-0.0016 |
-3.3531e-006 |
-3.3531e-006 |
0.5 |
0.0026 |
0.0026 |
5.4360e-006 |
5.4360e-006 |
0.6 |
0.0038 |
0.0016 |
-8.2015e-006 |
-3.3531e-006 |
0.7 |
-0.0053 |
-9.1897e-004 |
-1.1823e-005 |
-1.8342e-006 |
0.8 |
-0.0071 |
-4.0588e-004 |
-1.6559e-005 |
-7.9939e-007 |
0.9 |
-0.0092 |
-1.0109e-004 |
-2.2794e-005 |
-1.9754e-007 |
加hann hanning卷积窗提高振幅校正精度还可用於fft/apfft法中.
“全相位FFT在VC下的测试结果和分析”一文中提出频率分辨率对精度的影响,在你的测试中,对f=49.3 当频率分辨率=1时,峰值频谱线在round(49)+1=50的位置, 如图一(a). 当频率分辨率=8时,峰值频谱线在round(49/8)+1=7的位置, 如图一(b).当频率分辨率=16时, 峰值频谱线在round(49/16)+1=4的位置, 如图一(c). 频率分辨率8和16时频谱线靠近零轴, 这时余弦信号的镜像频率(图中频率轴[0 N]右端)的泄漏严重影响精度, 特别你用fft/apfft法,fft的泄漏影响大.所以你测的频率分辨率对精度的影响,主要是余弦信号的镜像频率泄漏的影响.
你将信号改为150赫测,结果就不一样, 将余弦信号改为单谱线的exp信号测,无噪时频率分辨率对校正精度没有影响.所以测试中频率分辨率一股选1, 频率分辨率>1相当密集频谱.
另外尽量使被测峰值运离[0 N/2]频谱两瑞靠近频谱中间.
GMT+8, 2024-5-15 11:41 , Processed in 0.032672 second(s), 16 queries , Gzip On.
Powered by Discuz! X3.4
Copyright © 2001-2021, Tencent Cloud.