第三章：信号的频域分析方法§3.1 频域分析的基本概念§3.2 傅氏变换§3.3 提高谱分析精度的一些方法§3.4 脉冲响应与传递函数§3.5 频谱的表示方法及频域特征提取方法§3.6 倒谱分析方法§3.7全息谱分析技术一.概述频谱：将动态信号的频率成分的幅值、相位、与频率的关系表达出来的图形。频谱有：离散谱——与周期性、准周期性信号对应；连续谱——与非周期信号及随机信号对应，用谱密度二.周期信号的频谱上式表明：周期分量可以分解为无限多个频率f为基频的整数倍的谐波分量之和。式中：0A0−−直流分量Akcos(2πkf0t+ϕk)－－谐波分量 为基频的整数倍Ak−−振幅ϕk－－相位f f10－－基频 0=T §3.1频域分析的基本概念一.概述频域分析是机械故障诊断中用得最广泛的信号处理方法之一。频谱分析：利用某种变换将复杂信号分解为简单信号的叠加。例：傅氏变换→把复杂信号分解为有限个或无限个频率的简谐分量；沃尔什变换→把复杂信号分解为不同的方波。二.周期信号的频谱周期信号：x(t)＝x(t+nT) 其中T——周期，n——整数。可展开成傅氏级数：∞x(t)=A0+∑Akcos(2πkf0t+ϕk)k=1(3－1)二.周期信号的频谱当周期信号只包含有限个谐波分量时，只有对应的Ak≠0，其余为0。通常采用谱图的形式来表示周期信号的谐波分量。最基本的谱图有：幅值谱——各谐波分量的振幅Ak与频率的图示关系；相位谱——各谐波分量的相位φk与频率的图示关系。1二.周期信号的频谱fff3f02f0X(t)A(f)f0φ(f)000周期信号的傅氏级数分解示意图二.周期信号的频谱二.周期信号的频谱傅氏级数也可以写成复指数函数形式。根据欧拉公式e±j2πft=cos2πft±jsin2πftcos2πft=1(e−j2πft+ej2πft2)sin2πft=j1(e−j2πft−ej2πft2)（3－1）可以表示成如下形式：x(t)=ej2πkf0t (k=0,±1,±2,...)n∑∞Ak=−∞(3－2)二.周期信号的频谱例：x(t)=5cos2Лf0t 只有一个谐波分量∴对应的幅值谱上只有一根谱线，相位谱上对应的相位为00x(t)=5cos2Лf0t＋3sin6Лf0t 有两个谐波分量∴对应的幅值谱上有两根谱线，相位谱上对应的相位为00、－900。二.周期信号的频谱各谐波分量叠加为复杂波形时，相位很重要。当各分量的Ak不变、φk改变时，会使合成后的信号波形变化很大。信号总能量为各谐波分量与直流分量能量之和，即x(t)的均方值为：x2rms=A20+122∑∞Akk=1二.周期信号的频谱其中，Ak是展开系数，其计算公式如下：TA1k=2T∫−Tx(t)e−j2πf0tdt（3－3）2式中：Ak为复数，由周期信号x(t)确定，它综合反映了n次谐波的幅值及初相信息2二.周期信号的频谱周期信号的频谱的特点：a.离散性b.谐波性c.收敛性周期信号的复频谱：其实频谱总是偶对称的；其虚频谱总是奇对称的。三.非周期信号的频谱振幅振幅A0t0Tt│X(f)││X(f)│10.50f0f01/T2/T3/Tfa.衰减振动波形及频谱b.半正弦脉冲波形及频谱非周期信号的波形及频谱三.非周期信号的频谱如果把f=ω2π代入式（3－4）、（3－5）可得：X(ω)=1∞2π∫+−∞x(t)e−jωtdt（3－6）x(t)=∫+∞jωt−∞X(ω)edω（3－7）X(f)=2πX(ω)（3－8）称X(f)为x(t)的连续频谱。一般X(f)是复函数。X(f)=X(f)ejφ(f)（3－9）三.非周期信号的频谱非周期信号（如瞬态信号、脉冲信号等）的频谱必须用连续谱表示。∵非周期信号可看作周期无穷大（T→∞）的周期信号，其基频趋向于0（f→0），因此，其谐波分量的间隔将无穷小，其频谱为连续谱。∴对非周期信号不能用幅值谱的概念，需要用谱密度的概念。三.非周期信号的频谱非周期信号x(t)的幅值谱密度和相位谱可以通过傅氏变换得到：X(f)=∫∞−∞x(t)e−j2πftdt（3－4）式中：j=−1虚数单位；尽管x(t)为实数，其傅氏变换X(f)一般为复数；其模X(f)为幅值谱密度；其幅角φ(f)=argX(f)为相位谱密度。X(f)的傅氏逆变换为：x(t)=∫∞−∞X(f)ej2πftdf（3－5）三.非周期信号的频谱式（3－9）代入（3－5）得：x(t)=∫∞ft+φ(f))−∞X(f)ej(2πdf（3－10）由： ej(2πft+ϕ(f))=cos(2πft+ϕ(f))+jsin(2πft+ϕ(f)) 取实部 得： x(t)＝∫∞－∞X(f)cos(2πft+ϕ(f))df )n与： x(t=A0＋∑Akcos(2πkf0t+ϕk) 相似k=1区别：│X(f)│代表幅值谱密度，是单位频带上的幅值，与信号幅值的量纲不一样。Ak代表幅值谱，量纲与信号幅值的量纲一样。3三.非周期信号的频谱X(t)与│X(f)│之间存在：∫∞∞−∞x2(t)dt=∫X(f)2−∞df (巴赛伐等式） (3－11)上式为总能量的频谱表达式，左边为X(t)在（－∞，＋∞）之间的总能量，右边│X(f)│2称为X(t)的能谱密度。三.非周期信号的频谱它的巴赛伐等式为:∫∞x2(t)dt=∫∞−∞X(f,T)2T−∞df (3－13)可得：lim1∞T→∞2T∫−∞x2T(t)dt=∫∞−∞lim1T→∞2TX(f,T)2df 上式左边为X(t)在（－∞，＋∞）上的平均功率，右边的被积式为平均功率谱密度，简称功率谱密度。记作：Sx(f)=lim1T→∞2TX(f,T)2(3－14)四.平稳随机信号的频谱1.自谱密度：S1[2x(f)=limT→∞TEX(f)](3－15)或：S)=lim12x(ωT→∞2πTE[X(ω)]自谱密度函数的性质：①Sx(f)是实的非负偶函数②Sx(f)是Rx(τ)的傅氏变换，Rx(τ)是Sx(f)的傅氏逆变换三.非周期信号的频谱∵许多时间函数（例如：正弦函数）的总能量无限，但其功率有限。∴考虑在（－∞，＋∞）上的平均功率：lim1Tx2T→∞2T∫−TT(t)dt x⎧⎪x(t) t≤TT(t)=⎨⎪⎩0 t>T 为x(t)的截尾函数其傅氏变换为：X(f,T)=∫TxT(t)e−j2πftdt=∫∞x2πft−T−∞T(t)e−jdt（3－12）四.平稳随机信号的频谱∵平稳随机信号不是周期信号∴其频谱应为连续谱又∵样本曲线的波形各不相同∴幅值谱没有意义∴平稳随机信号的频谱是指功率谱密度。四.平稳随机信号的频谱即：⎧j2⎪Sπfτx(f)=dτ ⎧⎪Sx(ω)=∞R(τ)⎨∫∞−∞Rx(τ)e−或⎨∫−∞xe−jωτdτ ⎪R∫∞⎩x(τ)=Sfτ⎪⎩R1∞x(τ)=−∞x(f)ej2πdf2π∫−∞Sx(ω)ejωτdω∵Sx(ω)、Rx(τ)均为实偶函数∴可用ejωτ的实部余弦函数来代替⎧⎪Sx(ω)=∫∞−∞Rx(τ)cosωτdτ ⎨⎪R(τ)=1∞(3－16)⎩x2π∫−∞Sx(ω)cosωτdω即：自谱密度与自相关函数是一个傅氏变换对。4四.平稳随机信号的频谱又有巴赛伐等式：信号时域中的总功率等于频域中的总功率。得：P=1∞x2(t)dt＝1∞X(ω21∞T∫−∞2πT∫−∞)dω =2π∫−∞Sx(ω)dω S2x(ω)=X(ω)/T (3－17)上式表明了功率谱与幅值谱的关系。四.平稳随机信号的频谱S(f)（a)窄带随机噪声的功率谱密度S(f)(c)白噪声的功率谱密度fS(f)fS(f)(b)宽带随机噪声的功率谱密度（d）正弦波加随机噪声的功率谱密度ff平稳随机信号的功率谱密度四.平稳随机信号的频谱Sxx(ω)的幅值为S2xy(ω)=Cxy(ω)+Q2xy(ω)相位：θω)xy=argtgQxy(Cxy[ω]互谱密度函数的标准化形式称作凝聚函数r2(ω)=S2xyxy(ω)/Sx(ω)Sy(ω)(3－20)0≤r2xy(ω)≤1四.平稳随机信号的频谱可见，要得到信号x(t)的功率谱有两条途径：相关图法——先求Rx(τ)，再求Sx(ω)周期图法——先求X(ω)，再求Sx(ω)四.平稳随机信号的频谱2.互谱密度S∞ωτxy(ω)=∫−∞Rxy(τ)e−jdτ互谱密度S)为复数。(3－18)xy(ωSxy(ω)=Cxy(ω)+jQxy(ω)（3－19）Cxy(ω)－－共谱密度Qxy(ω)－－重谱密度§3.2傅氏变换一.傅氏变换的基本性质X(f)=∞（3－21）若∫j2πft−∞x(t)e−dtx(t)=∫∞X(f)ej2πft（3－22）−∞dt则：称x(t)与X(f)为傅氏变换对。记作：x(t)←→X(f)X(f)为x(t)的傅氏变换；x(t)为X(f)的傅氏逆变换。5一.傅氏变换的基本性质1.线性叠加性若：x1(t)、x2(t)的傅氏变换为：X1(f)、X2(f)则：x1(t)＋x2(t)的傅氏变换为：X1(f)＋X2(f) ∫∞x−j2πft∞−j2πft∫∞−j2πft−∞(x1(t)+2(t))edt=∫−∞x1(t)edt+−∞x2(t)edt=X1(f)+X2(f)一般：C1x1(t)＋C2x2(t)←→C1X1(f)＋C2X2(f)C1、C2为常数。二.傅氏变换3.尺度变换如果：x(t)←→X(f) 令：t’=kt，其中k为大于0则：的常数x(kt)←→1⎛f⎞kX⎜⎝k⎟⎠∫∞x(kt)e−j2π2tdt=∫∞πfkt−∞x(t'−∞)e−j2d⎛⎜t'⎞1⎜⎛f⎞k⎟⎟=kX⎝⎠⎜⎝k⎟⎠即：时间尺度扩展或压缩k。k倍，相应的频率尺度压缩到1/k或扩展随时间尺度的压缩或扩展，频域的幅值也扩大或缩小k倍。所以：频域曲线下的面积保持不变。同理：1kx⎛⎜t⎞⎝k⎟⎠←→X(kf)一.傅氏变换的基本性质5. 频移定理x(t)ej2πf0t←→X(f−f0)如果X(f)的自变量移动一个常量f0，其傅氏逆变换x(t)要乘以ej2πf0t，产生调制现象，即x(t)与cos2лf0t相乘。x(t)tfx(t)cos2лf0t[X(f-f0)+X(f+f0)]/2tf0二.傅氏变换2对称性如果：x(t)←→X(f)则：x(±t)↔X(mf)也就是说：如果X(f)是信号x(t)的谱，则：x(±t)的谱是X(mf)。一.傅氏变换的基本性质4.时移定理如果：x(t)←→X(f) 则：x(t−t0)←→X(f)e−j2πft0令：s＝t－t0∫∞−∞x(t−tj2πft0)e−dt=∫∞−∞x(s)e−j2πf(s+t0)ds=∫∞−∞x(s)e−j2πfse−j2πft0ds=e−j2πft∞0∫j2πfs−∞x(s)e−ds即：时间位移引起相角=e−j2πftΦ(f)的变化，但不改变0X(f)傅氏变换的大小。一.傅氏变换的基本性质6. 卷积与乘积信号x1(t)、x2(t)的卷积记作：x1(t)*x2(t)，并有：x1(t)∗x2(t)=∫∞−∞x1(τ)x2(t−τ)dτ定义：如果：x1(t)←→X1(f)、x2(t)←→X2(f)则：x1(t)*x2(t)←→X1(f)X2(f)即：卷积的傅氏变换为各自傅氏变换的乘积。反之：x1(t) x2(t)←→X1(f)* X2(f)乘积的傅氏变换为各自傅氏变换的卷积。6二.离散傅氏变换(DFT)1.离散傅立叶变换的过程a. 对连续信号采样变成离散信号，为避免频率混淆，取fs≥2fmaxb.对采样值加窗截断，得到有限个采样数据，为减少频率泄漏必须使采样长度T0=NΔt足够长。c. 离散信号N个近似值对应X(f)的N个近似值，构成N个傅氏变换对。二.离散傅氏变换(DFT)∵采样间隔的大小不影响离散傅氏变换的实质；∴可以略去采样间隔，式3－23、3－24可以表示成：N−1X(n)=∑x(k)WnkN (n=0,1,2,...,N−1)（3－25）k=0N−1x(k)=1N∑X(n)W−nkN (k=0,1,2,...,N−1)n=0（3－26）式中：W=e−j2π/NN三.快速傅氏变换(FFT)因为DFT计算所需的时间与长度的平方成正比，所以，FFT算法的计算时间可缩短N/log2N。FFT与DFT算法的速度比为：N2/Nlog2N=N/log2N（3－27）序列长度N416642561024409616384速度比2410.732102.4341.31170.3二.离散傅氏变换(DFT)2. 离散傅氏变换的算法离散傅氏变换的表达式为N−1Χ(nΔf)=∑x(kΔt)e−j2πnk/Nn=0,1,LN−1k=0(3－23)Χ(nΔf)为复数幅值：Χ(nΔf)=R2(nΔf)+I2(nΔf)相角：ϕ(nΔf)=arctg⎢⎡I(nΔf)⎤⎣R(nΔf)⎥⎦逆变换：x(kΔt)=1∑N−1Χ(nΔf)ej2πnk/NNK=0,1, …,N-1n=0(3－24)三.快速傅氏变换(FFT)在用DFT方法计算N个X(nΔf)的值时，计算整个序列需作N2次复数乘法与N（N－1）加法运算，一次复数乘法相当于4次实数乘法，一次复数加法相当于2次实数加法，当N大时，DFT算法的计算量很大。三.快速傅氏变换(FFT)1. FFT算法的基本思想把长序列的DFT分解成几个短序列的DFT，并利用WknN的周期性和对称性减少DFT运算次数。最常用的是基2FFT( N=2M的FFT)FFT算法分：时域抽取法FFT (Decimation－In－Time FFT，简称DIT－FFT)；频域抽取法FFT (Decimation－In－Frequency FFT，简称DIF－FFT)。7三.快速傅氏变换(FFT)2. 时域抽取法基2FFT算法的过程(以N=8为例)先对原数据序列按奇偶逐步进行分解。原始序列：x0x1x2x3x4x5x6x71个长度为8的序列第一次分解: x0x2x4x6x1x3x5x72个长度为4的序列第二次分解: x0x4x2x6x1x5x3x74个长度为2的序列第三次分解: x0x4x2x6x1x5x3x78个长度为1的序列三.快速傅氏变换(FFT)又∵WN/2=e−j(2π/N)⋅N/2N=−1∴Wk+N/2N=WkN/2kN⋅WN=−WNX(k)=G(k)+WkNH(k) k=0,1,...,N/2−1（3－30）X(k+N/2)=G(k)−WkNH(k) k=0,1,...,N/（2−13－31）将两个半段X(k)和X(k+N/2)合成得到整个序列的X(k)。在合成时：偶序列DFT的变换G(k)不变，奇序列DFT 的变换H(k)要乘以权重函数WkN。二者合成时前半段的用加，后半段的用减。三.快速傅氏变换(FFT)X(0)X(1)W0X(2)X(3)WX(0)X(1)XW0(2)XW(3)N点DFT的第二次时域抽取分解图（N＝8）三.快速傅氏变换(FFT)由上述抽取方法及FFT的计算公式(3－25)可得：N/2−1X(k)=∑[x(2n)W2nk+x(2n+1)W(2n+1)kNN] k=0,1,...,N−1n=0（3－28）∵W2−2j(2π/N)N=e=e−j2π/(N/2)=W1N/2N/2−1X(k)=∑[x(2n)WnkN/+x(2n+1)WnkWk2N/2N] k=0,1,...,N−1n=0 =G(k)+WkNH(k)（3－29）式中：N∑/2−1N/2−1G(k)=x(2n)WnkN/2, H(k)=1)WnkN/2 k=0,1,...,N−1n=0∑x(2n+n=0G(k)、H(k)的周期为N/2即：G(k)=G(k+N/2)H(k)=H(k+N/2)三.快速傅氏变换(FFT)x(0)x1(0)X(0)x(1)x1(1)X(1)x(2)x1(2)x(3)xX(2)1(3)X(3)x(4)x2(0)x(5)xW0)2(1)X(4W2)X(5)x(6)x2(Wx(7)xX(6)2(3)WX(7)N点DFT的第一次时域抽取分解图（N＝8）三.快速傅氏变换(FFT)A(0)A(0)A(0)A(0)A(1)A(1)A(2)A(2)A(3)A(3)A(4)A(4)A(5)A(5)A(6)A(6)A(7)A(7)A(7)A(7)N点DIT－FFT运算流程（N＝8）8三.快速傅氏变换(FFT)N−1FFT逆变换的计算x(n)=1N∑X(k)W−knNn=0也可以按照上述方法进行。这时需将x(n)改为X(n)，并将系数1/N分解为m个1/2的连乘，每一级迭代乘上一个1/2因子。一.窗函数的选择对x(t)采一段样，相当于用矩形窗去截取原信号：当原始信号为非周期信号或者当原始信号为周期信号但截取长度不等于整周期时，就改变了原始信号的频率结构。谱图上会出现原频率以外的许多频率成分，即出现了频率泄漏现象。 整周期截断无泄漏非整周期截断有泄漏非整周期截断引起频率泄漏一.窗函数的选择2.窗函数的指标：①最大旁瓣值:用最大旁瓣峰值与主瓣峰值之比取对数表达20lg（A旁/A主）dB,此值为负数，越小越好。②旁瓣衰减率: 10个相邻旁瓣峰值的衰减比的对数表示，记作dB/10 oct,比值大衰减快，泄漏少。③主瓣宽：以下降3dB时的带宽表示，用3dB带宽乘以△f，△f表示分辨率。主瓣窄则可精确确定出峰值频率。§3.3提高谱分析精度的一些方法从傅氏变换的过程我们可以看出存在一些固有的缺陷，它影响了谱分析的精度，主要是：（1）用矩形窗去截断采样信号，引起了频率泄漏，使得X(f)产生误差；（2）因数据长度有限，对应点数有限，使频率分辨率受到限制；（3）平稳随机信号的傅氏变换受到样本曲线随机性的影响，误差较大。常用的提高频谱分析精度的方法有：窗函数、频率细化、平均等。一.窗函数的选择1. 常用的窗函数①Hanning窗（汉宁窗）w(t)={0.5−0.5cosϖt.........0≤t≤T0......t>0,t>T（3－32）②Hamming窗（海明窗）w(t)={0.54−0.85cosϖ,tKKt<0,t>TKK0≤t≤T0,（3－33）③指数窗w(t)={e−αt/T,tL<0L0≤t≤T0,LL,t>T（3－34)④矩形窗w(t)={1,0L,LL0≤t≤TLt<0,t>T(3－35)一.窗函数的选择④主瓣顶点最大误差：以％表示a 矩形窗泄漏严重，主瓣顶点误差最大，主瓣最窄，用于需精确确定出主瓣峰值频率时；b 汉宁窗泄漏最少，使用最多但主瓣较宽，主瓣误差也较大，需修正；c 指数窗无旁瓣，但主瓣明显加宽，主要用于脉冲或衰减信号，以提高信噪比。常用窗函数的对比窗函数最大旁瓣旁瓣衰减主瓣宽主瓣顶点最大误值（dB）率（dB/10oct）3dB带宽乘以△f差（％）矩形－13（21％）－200．89－36.34Hanning-13.47(2.7%)-601.44-15.12Hamming-43.19(0.7%)-201.3418.14指数－－宽9矩形窗函数矩形窗函数幅频曲线W(f)Tf-11T0T汉宁窗函数汉宁窗函数幅频曲线二.频率细化技术设：原频率分析范围为fmax，采样点数为N，fs＝2.56fmax，则：频率分辨率△f＝2.56fmax/N若要在f1～f2区间内将△f提高M倍，f2-f1≤fmax/M ，则采样长度应为MN点，具体步骤：二.频率细化技术⑷对Zt每隔M个抽选一个组成Wt，Wt=ZMt，t=0,1,2,…..N-1;⑸对Wt进行N点FFT得细化频谱，即在f2-f1范围内的分辨率为Δf΄=Δf/M，改变f2-f1，可重新按上述过程进行。二.频率细化技术细化谱的作用：①提高频率分辨率；②提高信噪比。频率细化（Zoom）来自照相技术的变焦距。从理论上讲提高频率分辨率只能增加采样长度，但这样会增加计算量，不易观察。这里简单介绍细化方法之一:复调制法二.频率细化技术⑴以采样频率fs采MN点，得xt，t＝0，1，2，……MN－1；⑵将序列xt乘以单位旋转矢量e-j2πf1t，得到新序列yt，对信号进行频移，将f1移到原点；⑶用截至频率为f2-f1的数字滤波器对yt进行数字滤波，得Zt，t=0，1，2，……MN－1；三.内插技术1. 谱线内插10三.内插技术2. 内插原理原信号是连续信号，而且从－∞到＋∞。采样：第一步：在时域里乘上一个脉冲，脉冲的间隔是采样间隔，如果采样频率为2KHz，时间间隔为0.005秒；第二步：截断，相当于原始信号乘上一个矩形函数；在时域里相乘相当于频域中的卷积，矩形函数在频域是：sinxx三.内插技术sin⎛)⎜Nω⎞ W(ω=⎝2⎟⎠exp⎡−j⎛N−1⎞ω⎤sin⎛⎜ω⎞⎢⎜2⎟⎠⎥⎝2⎟⎣⎝⎦⎠式中： ω－圆频率； N－采样长度。 ω=2πNk, k=1,2,L,N2−1或： W(k)=sin(πk)exp[−jkπ]≈Nsin(πk)exp[−jkπ]sin⎛⎜πk⎞πk⎝N⎟⎠对于矩形窗，其频谱主瓣中心恰为两相距单位距离的谱线的重心位置。相位exp(-jkπ) 正比于k，当变为k+1时，相位从0度变为180度。三.内插技术4. 谱峰幅值插值第i阶谱峰幅值为： A2πδAlefti=Nsin(δπ);或 A2π(1−δ)Arighti=Nsin(1−δ)π。三.内插技术∴采下来的是1倍频分量，实际上是sinxx内的两个点，现在已知这两个点及曲线方程，可以求出中间的极值点。三.内插技术3. 谱峰内插设：Aleft和Aright 分别为第i阶谱线的左右幅值；Fleft与Fright 分别为左右谱线的频率位置; N 为采样长度;ts为采样周期;则第i阶谱峰频率为： F⎛δ⎞⎛1−δ⎞i=Fleft+⎜⎝N•t⎟； 或 Fi=Fright+⎜⎟；s⎠⎝N•ts⎠式中： δ=⎛⎜Aright⎞⎜A⎟⎟； 窗为矩形窗；⎝left+Aright⎠ δ=⎛⎜2Aright−Aleft⎞⎜⎟⎝Aright+Aleft⎟； 窗为海宁窗。⎠三.内插技术5. 相位插值初始时间相位：θi=phase⎡⎣Aleft⎤⎦−Aiδ+0.5πθi=phase⎡⎣Aright⎤⎦−Ai(δ−1)+0.5π在一个轴承断面上的X和Y分量相同谐波之间的相对相位是准确的，不需要任何修正。11三.内插技术5. 相位插值信号中第i阶谐波相对于基频分量的相位为：θ−θphase⎡⎣Ai1i1=left⎤⎦−iphase⎡⎣Aleft⎤⎦+α式中α的值可如下选择：iδ<1， α=0.5(i−1)π；1