时域和频域特征提取Matlab编程实例

1、信号处理课程设计第一章 绪论1.1 概述机械信号是指机械系统在运行过程中各种随时间变化的动态信息，经各种测试仪器拾取并记录和存储下来的数据或图像。机械设备是工业生产的基础，而机械信号处理与分析技术则是工业发展的一个重要基础技术。随着各行各业的快速发展和各种各样的应用需求，信号分析和处理技术在信号处理速度、分辨能力、功能范围以及特殊处理等方面将会不断进步，新的处理激素将会不断涌现。当前信号处理的发展主要表现在：1.新技术、新方法的出现；2.实时能力的进一步提高；3.高分辨率频谱分析方法的研究三方面。信号处理的发展与应用是相辅相成的，工业方面应用的需求是信号处理发展的动力，而信号处理的发展反过来又

2、拓展了它的应用领域。机械信号的分析与处理方法从早期模拟系统向着数字化方向发展。在几乎所有的机械工程领域中，它一直是一个重要的研究课题。机械信号分析与处理技术正在不断发展，它已有可能帮助从事故障诊断和监测的专业技术人员从机器运行记录中提取和归纳机器运行的基本规律，并且充分利用当前的运行状态和对未来条件的了解与研究，综合分析和处理各种干扰因素可能造成的影响，预测机器在未来运行期间的状态和动态特性，为发展预知维修制度、延长大修期及科学地制定设备的更新和维护计划提供依据，从而更为有效地保证机器的稳定可靠运行，提高大型关键设备的利用率和效率。机械信号处理是通过对测量信号进行某种加工变换，削弱机械信号中的

3、无用的冗余信号，滤除混杂的噪声干扰，或者将信号变成便于识别的形式以便提取它的特征值等。机械信号处理的基本流程图如图1.1所示。 测量信号信号处理 特征提取获得设备结构特征或运行状态机械设备图1.1 机械信号处理的基本流程 本文主要就第三、第四步骤展开讨论。第2章 机械信号的时域处理及其分析方法2.1 时域统计特征参数处理通过时域波形可以得到的一些特征参数，它们常用于对机械进行快速评价和简易诊断。2.1.1 有量纲的幅值参数有量纲的幅值参数包括方根幅值、平均幅值、均方幅值和峰值等。若随机过程x(t)符合平稳、各态历经条件且均值为零，设x为幅值，p(x)为概率密度函数，有量纲型幅值参数可定义为xd

4、=式中：xr为方根均值，为均值，为均方值，为峰值。由于有量纲型幅值参数来描述机械状态，不但与及其的状态有关，而且与机器的运动参数（如转速、载荷等）有关，因此直接用它们评价不同工况的机械无法得出统一的结论。2.1.2 无量纲型参数无量纲型参数具有对机械工况变化不敏感的特点，这就意味这，理论上它们与机械的运动条件无关，它们只依赖于概率密度函数p(x)的形状，所以无量纲型参数是一种较好的评价参数。一般它可定义为，由此公式，可得到如下的一些指标波形指标l=2,m=1 K=峰值指标l,m=2 C=脉冲指标l,m=1 I=裕度指标l,m=1/2 L=峭度指标 K=式中为信号标准差=2.2 相关分析方法以及

5、应用所谓相关，就是指变量之间的线性关系，它是一个非常重要的概念。对于确定性信号，两个变量之间可以用函数关系来描述，两者一一对应并为确定的数值。而两个随即变量之间不具有确定的关系。但是，如果这两个变量之间存在着某种不确定但却有着表征其特性的近似关系，这两个变量之间会有一定的线性关系。这时，对于一个随机机械信号，可以采用相关性函数来描述其在不同时间的幅值变化相关程度。2.2.1 自相关函数的概念和性质 x(t)是各态历经随机过程的一个样本函数，x(t+t)是x(t)时移t后的样本(图2.6)，把相关系数rx(t)x(t+t)简写为rx(t)，那么就有: 图2.6 波形图若用Rx(t)表示自相关函数

6、，其定义为: 信号的性质不同，自相关函数有不同的表达形式。如对周期信号（功率信号）：非周期信号（能量信号）： 图2.7给出了自相关函数具有的性质。正弦函数的自相关函数是一个余弦函数，在=0时具有最大值。它保留了幅值信息和频率信息,但丢失了原正弦函数中的初始相位信息。2.3 Matlab编程实验结果2.3.1 构造加噪周期信号，时域特征分析，自相关函数特性的验证，（程序1）图2.8 噪声-自相关.jpg 如图所示：自相关函数消除了大量的噪声，周期成分变得非常明显。原始信号的时域处理结果：平均值：0.0184极小值：-2.8138极大值：2.8557标准差：1.0103 方差： 1.0207峰峰值

7、：5.6695第3章 机械信号的频域处理方法及其应用信号处理中，傅立叶变换把一个随机信号解析成不同频率的正弦波，使信号的频域分析称为可能。由于计算机技术的发展，在微机上直接使用离散傅立叶变换变得非常方便，这使得频域分析称为常用的处理方法。常用的频域分析方法包括自谱、功率谱、倒谱等。3.1 频谱的分析方法DFT和FFT3.1.1 离散傅立叶变换DFT傅立叶变换及其逆变换都不适合用数字计算机计算。要进行数字计算和处理，必须将连续信号离散化，无限数据有限化。这种对有限个离散数据的傅立叶变换，称为有限离散傅立叶变换，简称DFT（Discrete Fourier Trasform）。3.1.2 快速傅立

8、叶变换FFT1965年J.W.Cooley 和J.W.Tukey研究一种DFT的快速算法，称为快速傅立叶变换，简称FFT(FastFourier Transform)。FFT的迅速发展，使数字频谱分析取得了突破性的进展。根据FFT快速变换的指导思想，就可以编制FFT的计算程序。时间序列从时域到频域要用FFT变换，从频域到时域要用逆变换IFFT，FFT和IFFT的公式可以统一。3.1.3 功率谱密度函数的物理意义Sx(f)和Sxy(f)是随机信号的频域描述函数。Sx(f)表示信号的功率密度沿频率轴的分布，故又称Sx(f)为功率谱密度函数。3.2 功率谱方法以及应用功率谱的定义式为若X（）=DFT

9、x(m)，x(n)为N点序列。则X（） =DFTx(-m)从而有 DFTR(M)= DFTx(m) DFTx(-m)即 （）= X（）X（）=| X（）|2综上所述，先用FFT求出随机离散序列的DFT，再计算幅频特性的平方，再除以N，即得到该随机信号的功率谱估计。3.3 倒频谱分析方法倒频谱实际上是频域信号取对数的傅立叶变换再处理，或称为“频域信号的傅立叶变换再变换”。对功率谱密度函数取对数的目的是使再变换以后信号的能量更加集中。倒频谱可以分析复杂频谱上的周期成分，分离和提取在密集泛频信号中的成分。对于具有同族谐频和异族谐频等复杂信号的分析，效果很好。倒频谱用于对语音分析中的语言音调的测定和检

10、测、机械振动谱图中的谐波分量作故障检测和诊断以及排除回波等方面是很有效的。3.3.1 倒频谱的数学描述倒频谱函数CF(q)(power cepstrum)其数学表达式为：CF(q)又叫功率倒频谱，或叫对数功率谱的功率谱。工程上常用的是式(2.67)的开方形式，即：C0(q)称为幅值倒频谱，有时简称倒频谱。倒频谱自变量q的物理意义为了使其定义更加明确，还可以定义：即倒频谱定义为信号的双边功率谱对数加权，再取其傅里叶逆变换，联系一下信号的自相关函数：看出，这种定义方法与自相关函数很相近，变量q与在量纲上完全相同。 为了反映出相位信息，分离后能恢复原信号，又提出一种复倒频谱的运算方法。若信号x(t)

11、的傅里叶变换为X(f): x(t)的倒频谱记为：显而易见，它保留了相位的信息。 倒频谱与相关函数不同的只差对数加权，目的是使再变换以后的信号能量集中，扩大动态分析的频谱范围和提高再变换的精度。还可以解卷积(褶积)成分，易于对原信号的分离和识别。3.4 细化谱分析方法细化谱分析法是增加频谱中某些部分分辨能力的方法，即“局部放大”的方法。所谓细化分析室只对固定某窄带部分进行放大，像照相机将照片的个别部分放大一样，使其动态范围和分辨率都提高。细化的分析过程中，首先像通常的FFT做法那样，选用采样频率fs=1/h进行采样，可得到N点离散序列xn.假设我们感兴趣的谱中心频率为fk的一个窄带f,然后用一个

12、复正弦序列（单位旋转矢量）exp-j2fknh乘以xn的yn新的N点离散序列。根据频移定理，即将频率原点有效地移至频率fk（即复调制）。fk成为新的频率坐标原点。正、负采样频率fs也同样移动了一个量fk。低通滤波后得到gm序列所保留下来的窄频带，若滤波后的总带宽小于采样频率的1/D倍，就有可能把采样频率降低到1/D，而不会再新的乃奎斯特频率附近产生混叠。然后再重新采样，用fs2= fs/D的频率来采样，即降低了采样频率。由采样定理可知，降低采样频率而又保持同样的采样点数N时，就相当于总的时间窗增长D倍，那么，频率分辨率也提高了D倍。所以，对经过重新采样后获得的新的离散序列rm进行复数FFT计算

13、，即可得到细化后的谱线，这些谱线就代表中心频率为fk的一窄带f间的细化谱。3.5 Matlab编程实验结果3.5.1 产生一组由60HZ和150HZ的正弦信号和随机噪声组成的信号，观察其时域波形与频谱。（程序2） 图3-1 原始信号的时域波形图 图3-2 原始信号的频谱图图3-1看不出信号 的周期成分；图3-2可以清除看到，在频率为60HZ和150HZ处有两个尖峰，即为信号的两个频率分量。3.5.2 功率谱估计（周期图法）：1.利用上图的带噪原始信号的傅里叶变换后结果幅值，将幅值平方，即可得功率谱的估计值（Welch法） 图3-3 采样点数为1024时的估计功率谱 图3-4 采样点数为256时

14、的估计功率谱由图3-3与3-4可看出：2.为提高周期图的平滑性，将信号分段估计并进行平均来减少功率谱估计的协方差，得到平均周期图。 图3-5 三段平均的估计功率谱图3-6 六段平均的估计功率谱由图3-5与3-6看出：分段平均法提高了功率谱图的平滑性，分段数越多，平滑效果越好，信号细节更易丢失。3.对数据分段加非矩形创形成修正的功率谱估计法：图3-7 加汉宁窗的估计功率谱由于窗在其边沿为零，这减少了分段对混叠的依赖效果。用合适的窗函数，采用分段长度一半的混叠率能极大地降低估计的协方差。3.5.3 倒频谱分析：图3-8 实倒谱图3-9 复倒谱 正弦信号，其第一个功率谱变换为一脉冲，经滤波后进入第二

15、次功率谱变换，其输出为幅度很低的三角波输出，因而检测不到其存在。3.5.4细化谱分析：图3-10 原始信号FFT图3-11 ZOOM-FFT程序清单程序1.构造加噪周期信号，时域特征分析，自相关函数特性的验证fs=1000;t=0:1/fs:(1-1/fs);maxlag=100;x=randn(1,fs);c,maxlags=xcorr(x,maxlag); %白噪声的自相关性z=cos(2*pi*20*t)+0.7*randn(1,1000);%加白噪声m= mean(z); disp (m); %计算平均值mi = min(z); disp (mi); %极小值mx = max(z);

16、disp (mx); %极大值st = std(z); disp (st); %标准差fc = st.2； %方差figure(1)subplot(2,2,1) %2*2第一张图plot(t,x) %图片区域大小xlabel(t);ylabel(x(t);title(白噪声); %加标题subplot(2,2,2)plot(maxlags/fs,c)xlabel(t);ylabel(r(t);title(白噪声自相关);c,lags=xcorr(z,maxlag); %带白噪声的余弦信号自相关subplot(2,2,3)plot(t,z)xlabel(t);ylabel(z(t);title(

17、原始信号);subplot(2,2,4)plot(maxlags/fs,c)xlabel(t);ylabel(rz(t); title(自相关); 程序2: 产生一组由60HZ和150HZ的正弦信号和随机噪声组成的信号，对其做频谱分析、倒谱分析以及几种种功率谱估计方法的比较。%1.(频谱分析)产生一组由60HZ和150HZ的正弦信号和随机噪声组成的信号，观察其时域波形与频谱。fs=1000;N=1024;t=(0:N-1)/fs;f1=60;f2=150;s1=sin(2*pi*f1*t)+sin(2*pi*f2*t);s2=2*randn(size(t);x=s1+s2;figure(1)s

18、ubplot(2,1,1)plot(t,x)X=abs(fft(x);f=(0:N/2-1)*fs/N;subplot(2,1,2)plot(f(1:N/2),X(1:N/2)%2.倒谱分析D=rceps(x); %实倒谱figure(2)subplot(2,1,1)plot(t,D) E=cceps(x); %复倒谱subplot(2,1,2)plot(t,E)%3.1功率谱估计（Welch法）Pxx=abs(fft(x,N).2/N; %采样点数为1024Pxx_short=abs(fft(x,256).2/256; %采样点数为256figure(3)subplot(2,1,1)plot

19、(0:N-1)/N*fs,10*log10(Pxx)subplot(2,1,2)plot(0:255)/256*fs,10*log10(Pxx_short)*10)%3.2将信号分段估计并进行平均来减少功率谱估计的协方差，得到平均周期图。Pxx=(abs(fft(x(1:256).2+abs(fft(x(257:512).2+abs(fft(x(513:768).2)/256/3;figure(4)plot(0:255)/256*fs,10*log10(Pxx) %3.3将信号分为六段作功率谱估计再平均。Pxx=(abs(fft(x(1:128).2+abs(fft(x(129:256).2+

20、abs(fft(x(257:384).2+abs(fft(x(385:512).2+abs(fft(x(513:640).2+abs(fft(x(641:768).2)/256/6;figure(4)plot(0:127)/128*fs,10*log10(Pxx)%3.4对数据分段家非矩形创形成修正的周期突法。窗在其边沿为零，这减少了分段对混叠的依赖效果。用合适的窗函数（如海明窗，汉宁窗），采用分段长度一半的混叠率能%极大地降低估计的协方差。汉宁法：w=hanning(256);Pxx=(abs(fft(w.*x(1:256).2+abs(fft(w.*x(129:384).2+abs(fft

21、(w.*x(257:512).2+abs(fft(w.*x(385:640).2+abs(fft(w.*x(513:768).2+abs(fft(w.*x(641:896).2)/(norm(w)2*6);figure(5)plot(0:255)/256*fs,10*log10(Pxx)程序3：%ZOOM-FFTfs=200;N=1024;n=0:N-1;t=n/fs;f=(0:N-1)*fs/N;f1=7;f2=7.2;f3=8;s1=sin(2*pi*t*f1);s2=sin(2*pi*t*f2);s3=sin(2*pi*t*f3);x=s1+s2+s3;load zoomfftdata;

22、fi=6;%最小细化截止频率np=10;%放大倍数nfft=512;%fft长度nt=length(x);fa=fi+0.5*fs/np;%最大细化截止频率nf=2nextpow2(nt);%？na=round(0.5*nf/np+1);%频移n=0:nt-1;b=n*pi*(fi+fa)/fs;%确定旋转因子y=x.*exp(-i*b);b=fft(y,nf);%fft变换a(1:na)=b(1:na);%正频率带通内的元素赋值a(nf-na+1:nf)=b(nf-na+1:nf);b=ifft(a,nf);%负频率带通内的元素赋值c=b(1:np:nt);%重采样a=fft(c,nfft)

23、*2/nfft;%进行ZOOM-FFT(nfft啥玩意儿？)%变换结果重新排序:y2=zeros(1,nfft/2);y2(1:nfft/4)=a(nfft-nfft/4+1:nfft);y2(nfft/4+1:nfft/2)=a(1:nfft/4);n=0:(nfft/2-1);f2=fi+n*2*(fa-fi)/nfft;%FFT变换y1=fft(x,nfft)*2/nfft;f1=n*fs/nfft;ni=round(fi*nfft/fs+1);na=round(fa*nfft/fs+1);%输出波形subplot(2,1,1)plot(t,x);subplot(2,1,2)nn=ni:na;plot(f1(nn),abs(y1(n