离散频谱校正理论和技术,不知道大家对这个名词熟不熟悉.近来在论坛上看到一些帖子讨论为何经FFT得到的幅值、频率和相位不准的。 其实前面我也发过一篇介绍离散频谱校正的综述性的文章,可能大家都忙,没时间去看,呵呵,这里我就我的理解,把离散频谱分析的误差来源和校正方法做个简单的介绍。另外,我想借这个机会了解下,在实际的应用过程中,非常精确的提取这些参数有没有必要,大家平常也没有碰到由于误差而带来的困扰,例如不能对故障做出判断,甚至出现错判的例子。大家平常用离散频谱校正用的多不多? 请大家讨论下,让我对我这个方向的应用也有个了解,谢谢! 离散频谱分析的误差产生的原因主要来自两方面,一方面是由于时域加窗截断产生的频域连续化,另一方面是由于计算机只能对有限的离散的频率进行计算,也即是频域离散化的结果。其中,加窗截断的影响使一个无穷长单频率信号在频域对应的一根谱线,变成一个连续谱,以加矩形窗为例,则是变成一个sinc型函数的形状,其峰值对应的频率即为单频信号的频率。但是由于频域的离散化,我们用FFT计算的频率一般都不会刚好会落在峰值处,这就是我们平时常说的泄露,这时我们就只能把计算得到的峰值谱线对应的频率做为估计的频率,如果以频率分辨率fs/N做归一 (即把频率分辨率看成1)的话,这个估计的频率的最大绝对值误差就是0.5,而幅值误差则依赖于加的窗的类型,由于矩形窗主瓣宽度为2,频谱开状较尖,幅值误差也就大。至于相位的最大误差则会相应的达到正负90度,已经完全不能用了。 离散频谱校正就是针对这种误差提出的各种校正出实际的频率、幅值和相位的一门理论和技术。国内现在比较常用的方法有比值(插值)法、能量重心法、FFT+FT法和相位差法,都有其各自的特点和优缺点。这里我给出一个比值校正法的程序供大家一起研究下。 当然,对于多频率成分的信号来说,离散频谱分析的另一个误差是来自于频率之间的相互干涉,这也是由于泄露所引起的,这个误差则主要靠加窗抑制旁瓣和减小频率分辨率、拉大频率间的距离(可通过ZFFT实现)来尽量减小。 %SpectrumCorrect_Test.m close all; clear all; clc; fs=1024; N=1024; t=(0:N-1)/fs; x=4*cos(2*pi*80*t+30*pi/180)+3*cos(2*pi*150.232*t+80*pi/180)+1*cos(2*pi*253.5453*t+240*pi/180); xf=fft(x); xf=xf(1:N/2)/N*2; XfCorrect=SpectrumCorrect(xf,3,1); XfCorrect(:,1)=XfCorrect(:,1)*fs/N; XfCorrect w=hann(N,'periodic'); xfw=fft(x.*w'); xfw=xfw(1:N/2)/N*4; XfCorrectW=SpectrumCorrect(xfw,3,2); XfCorrectW(:,1)=XfCorrectW(:,1)*fs/N; XfCorrectW %离散频谱比值校正法 by yangzj 2007.4.28 % %xf为FFT后的复数谱 %CorrectNum为校正的谱线条数,即校正最大的CorrectNum条 %WindowType为加窗类型,1为矩形窗,2为Hanning窗 % %SpectrumCorrect.m function XfCorrect=SpectrumCorrect(xf,CorrectNum,WindowType) XfCorrect=zeros(CorrectNum,3); for i=1:CorrectNum A=abs(xf); [Amax,index]=max(A); phmax=angle(xf(index)); %比值法 %加矩形窗 if (WindowType==1) indsecL=A(index-1)>A(index+1); df=indsecL.*A(index-1)./(Amax+A(index-1))-(1-indsecL).*A(index+1)./(Amax+A(index+1)); XfCorrect(i,1)=index-1-df; XfCorrect(i,2)=Amax/sinc(df); XfCorrect(i,3)=(phmax+pi*df)*180/pi; xf(index-2:index+2)=zeros(1,5); end %比值法 %加Hanning窗 if (WindowType==2) indsecL=A(index-1)>A(index+1); df=indsecL.*(2*A(index-1)-Amax)./(Amax+A(index-1))-(1-indsecL).*(2*A(index+1)-Amax)./(Amax+A(index+1)); XfCorrect(i,1)=index-1-df; XfCorrect(i,2)=(1-df^2)*Amax/sinc(df); XfCorrect(i,3)=(phmax+pi*df)*180/pi; xf(index-4:index+4)=zeros(1,9); end XfCorrect(i,3)=mod(XfCorrect(i,3),360); XfCorrect(i,3)=XfCorrect(i,3)-(XfCorrect(i,3)>180)*360; end 运行结果 XfCorrect = 80.0014 4.0016 29.8261 150.2333 2.9981 79.7127 253.5397 0.9996 -118.7272 XfCorrectW = 80.0000 4.0000 30.0000 150.2320 3.0000 80.0000 253.5453 1.0000 -120.0002 [ 本帖最后由 yangzj 于 2007-4-28 15:43 编辑 ] |
原帖由 yangzj 于 2007-4-30 12:37 发表
呵呵,看来是没有市场哦
原帖由 w89986581 于 2007-5-6 15:36 发表
我尝试用全相法修正离散频谱,然后估计相位,没有成功过:(
原帖由 zhwang554 于 2007-5-7 12:50 发表
用全相位测相位的程序在本论坛的"请教全相位谱分析问题"和"如何准确确定信号中强线谱的相位?"有介绍
非常精确的提取信号的频率、幅值和相位参数的意义是很大的,用FFT非常精确的提取信号的频率、幅值和相位参数 ...
原帖由 zhwang554 于 2007-5-7 15:12 发表
由於数学公式贴不上,简单说一下.
全相位FFT的泄漏是原FFT泄漏的平方,泄漏分贝减少一倍,如FFT泄漏-20db,全相位
FFT泄漏为-40db.你用apFFT的程序将apFFT的振幅谱和原FFT的振幅谱画出来,一比就知道了.
复 ...
原帖由 shenyongjun 于 2007-5-7 17:02 发表
丁康老师有一篇关于频谱校正的综述论文发表在几年前的《振动工程学报》,最近又在看丁老师《齿轮箱故障诊断实用技术》(书名记不清了)。不知版主是否了解丁老师的工作?
原帖由 zhwang554 于 2007-5-7 17:34 发表
你这个问题提得好由克拉美-罗限而言,N增大克拉美-罗限要求更高.
N阶apFFT需要2N-1个数据,其泄漏是N阶FFT的平方,
用这2N-1个数据作FFT,其泄漏仍差於N阶apFFT,见附图.
而且2N-1个数据作FFT的计算量大了1倍.
...
GMT+8, 2024-11-26 06:01 , Processed in 0.077639 second(s), 25 queries , Gzip On.
Powered by Discuz! X3.4
Copyright © 2001-2021, Tencent Cloud.