时域和频域特征提取Matlab编程实例.docx
《时域和频域特征提取Matlab编程实例.docx》由会员分享，可在线阅读，更多相关《时域和频域特征提取Matlab编程实例.docx（20页珍藏版）》请在冰点文库上搜索。
时域和频域特征提取Matlab编程实例
第一章绪论
1.1概述
机械信号是指机械系统在运行过程中各种随时间变化的动态信息，经各种测试仪器拾取并记录和存储下来的数据或图像。
机械设备是工业生产的基础，而机械信号处理与分析技术则是工业发展的一个重要基础技术。
随着各行各业的快速发展和各种各样的应用需求，信号分析和处理技术在信号处理速度、分辨能力、功能范围以及特殊处理等方面将会不断进步，新的处理激素将会不断涌现。
当前信号处理的发展主要表现在：
1.新技术、新方法的出现；2.实时能力的进一步提高；3.高分辨率频谱分析方法的研究三方面。
信号处理的发展与应用是相辅相成的，工业方面应用的需求是信号处理发展的动力，而信号处理的发展反过来又拓展了它的应用领域。
机械信号的分析与处理方法从早期模拟系统向着数字化方向发展。
在几乎所有的机械工程领域中，它一直是一个重要的研究课题。
机械信号分析与处理技术正在不断发展，它已有可能帮助从事故障诊断和监测的专业技术人员从机器运行记录中提取和归纳机器运行的基本规律，并且充分利用当前的运行状态和对未来条件的了解与研究，综合分析和处理各种干扰因素可能造成的影响，预测机器在未来运行期间的状态和动态特性，为发展预知维修制度、延长大修期及科学地制定设备的更新和维护计划提供依据，从而更为有效地保证机器的稳定可靠运行，提高大型关键设备的利用率和效率。
机械信号处理是通过对测量信号进行某种加工变换，削弱机械信号中的无用的冗余信号，滤除混杂的噪声干扰，或者将信号变成便于识别的形式以便提取它的特征值等。
机械信号处理的基本流程图如图1.1所示。
图1.1机械信号处理的基本流程
本文主要就第三、第四步骤展开讨论。
第2章机械信号的时域处理及其分析方法
2.1时域统计特征参数处理
通过时域波形可以得到的一些特征参数，它们常用于对机械进行快速评价和简易诊断。
2.1.1有量纲的幅值参数
有量纲的幅值参数包括方根幅值、平均幅值、均方幅值和峰值等。
若随机过程x（t）符合平稳、各态历经条件且均值为零，设x为幅值，p（x）为概率密度函数，有量纲型幅值参数可定义为
xd=
=
式中：
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相关分析方法以及应用
所谓相关，就是指变量之间的线性关系，它是一个非常重要的概念。
对于确定性信号，两个变量之间可以用函数关系来描述，两者一一对应并为确定的数值。
而两个随即变量之间不具有确定的关系。
但是，如果这两个变量之间存在着某种不确定但却有着表征其特性的近似关系，这两个变量之间会有一定的线性关系。
这时，对于一个随机机械信号，可以采用相关性函数来描述其在不同时间的幅值变化相关程度。
2.2.1自相关函数的概念和性质
x（t）是各态历经随机过程的一个样本函数，x（t+）是x（t）时移后的样本（图2.6），把相关系数x（t）x（t+）简写为x（），那么就有:
图2.6波形图
若用Rx（）表示自相关函数，其定义为:
信号的性质不同，自相关函数有不同的表达形式。
如对周期信号（功率信号）：
非周期信号（能量信号）：
图2.7给出了自相关函数具有的性质。
正弦函数的自相关函数是一个余弦函数，在τ=0时具有最大值。
它保留了幅值信息和频率信息,但丢失了原正弦函数中的初始相位信息。
2.3Matlab编程实验结果
2.3.1构造加噪周期信号，时域特征分析，自相关函数特性的验证，（程序1）
图2.8噪声--自相关.jpg
如图所示：
自相关函数消除了大量的噪声，周期成分变得非常明显。
原始信号的时域处理结果：
平均值：
0.0184
极小值：
-2.8138
极大值：
2.8557
标准差：
1.0103
方差：
1.0207
峰峰值：
5.6695
第3章机械信号的频域处理方法及其应用
信号处理中，傅立叶变换把一个随机信号解析成不同频率的正弦波，使信号的频域分析称为可能。
由于计算机技术的发展，在微机上直接使用离散傅立叶变换变得非常方便，这使得频域分析称为常用的处理方法。
常用的频域分析方法包括自谱、功率谱、倒谱等。
3.1频谱的分析方法
DFT和FFT
3.1.1离散傅立叶变换DFT
傅立叶变换及其逆变换都不适合用数字计算机计算。
要进行数字计算和处理，必须将连续信号离散化，无限数据有限化。
这种对有限个离散数据的傅立叶变换，称为有限离散傅立叶变换，简称DFT（DiscreteFourierTrasform）。
3.1.2快速傅立叶变换FFT
1965年J.W.Cooley和J.W.Tukey研究一种DFT的快速算法，称为快速傅立叶变换，简称FFT（FastFourierTransform）。
FFT的迅速发展，使数字频谱分析取得了突破性的进展。
根据FFT快速变换的指导思想，就可以编制FFT的计算程序。
时间序列从时域到频域要用FFT变换，从频域到时域要用逆变换IFFT，FFT和IFFT的公式可以统一。
3.1.3功率谱密度函数的物理意义
Sx（f）和Sxy（f）是随机信号的频域描述函数。
Sx（f）表示信号的功率密度沿频率轴的分布，故又称Sx（f）为功率谱密度函数。
3.2功率谱方法以及应用
功率谱的定义式为
若X（Ω）=DFT[x（m）]，x（n）为N点序列。
则
X
（Ω）=DFT[x
（-m）]
从而有DFT[R（M）]=
DFT[x（m）]DFT[x
（-m）]
即
（Ω）=
X（Ω）X
（Ω）=
|X（Ω）|^2
综上所述，先用FFT求出随机离散序列的DFT，再计算幅频特性的平方，再除以N，即得到该随机信号的功率谱估计。
3.3倒频谱分析方法
倒频谱实际上是频域信号取对数的傅立叶变换再处理，或称为“频域信号的傅立叶变换再变换”。
对功率谱密度函数取对数的目的是使再变换以后信号的能量更加集中。
倒频谱可以分析复杂频谱上的周期成分，分离和提取在密集泛频信号中的成分。
对于具有同族谐频和异族谐频等复杂信号的分析，效果很好。
倒频谱用于对语音分析中的语言音调的测定和检测、机械振动谱图中的谐波分量作故障检测和诊断以及排除回波等方面是很有效的。
3.3.1倒频谱的数学描述
倒频谱函数CF（q）（powercepstrum）其数学表达式为：
CF（q）又叫功率倒频谱，或叫对数功率谱的功率谱。
工程上常用的是式（2.67）的开方形式，即：
C0（q）称为幅值倒频谱，有时简称倒频谱。
倒频谱自变量q的物理意义
为了使其定义更加明确，还可以定义：
即倒频谱定义为信号的双边功率谱对数加权，再取其傅里叶逆变换，联系一下信号的自相关函数：
看出，这种定义方法与自相关函数很相近，变量q与τ在量纲上完全相同。
为了反映出相位信息，分离后能恢复原信号，又提出一种复倒频谱的运算方法。
若信号x（t）的傅里叶变换为X（f）:
x（t）的倒频谱记为：
显而易见，它保留了相位的信息。
倒频谱与相关函数不同的只差对数加权，目的是使再变换以后的信号能量集中，扩大动态分析的频谱范围和提高再变换的精度。
还可以解卷积（褶积）成分，易于对原信号的分离和识别。
3.4细化谱分析方法
细化谱分析法是增加频谱中某些部分分辨能力的方法，即“局部放大”的方法。
所谓细化分析室只对固定某窄带部分进行放大，像照相机将照片的个别部分放大一样，使其动态范围和分辨率都提高。
细化的分析过程中，首先像通常的FFT做法那样，选用采样频率fs=1/h进行采样，可得到N点离散序列{xn}.假设我们感兴趣的谱中心频率为fk的一个窄带
f,然后用一个复正弦序列（单位旋转矢量）exp[-j2
fknh]乘以{xn}的{yn}新的N点离散序列。
根据频移定理，即将频率原点有效地移至频率fk（即复调制）。
fk成为新的频率坐标原点。
正、负采样频率
fs也同样移动了一个量fk。
低通滤波后得到{gm}序列所保留下来的窄频带，若滤波后的总带宽小于采样频率的1/D倍，就有可能把采样频率降低到1/D，而不会再新的乃奎斯特频率附近产生混叠。
然后再重新采样，用fs2=fs/D的频率来采样，即降低了采样频率。
由采样定理可知，降低采样频率而又保持同样的采样点数N时，就相当于总的时间窗增长D倍，那么，频率分辨率也提高了D倍。
所以，对经过重新采样后获得的新的离散序列{rm}进行复数FFT计算，即可得到细化后的谱线，这些谱线就代表中心频率为fk的一窄带
f间的细化谱。
3.5Matlab编程实验结果
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时的估计功率谱
由图3-3与3-4可看出：
2.为提高周期图的平滑性，将信号分段估计并进行平均来减少功率谱估计的协方差，得到平均周期图。
图3-5三段平均的估计功率谱
图3-6六段平均的估计功率谱
由图3-5与3-6看出：
分段平均法提高了功率谱图的平滑性，分段数越多，平滑效果越好，信号细节更易丢失。
3.对数据分段加非矩形创形成修正的功率谱估计法：
图3-7加汉宁窗的估计功率谱
由于窗在其边沿为零，这减少了分段对混叠的依赖效果。
用合适的窗函数，采用分段长度一半的混叠率能极大地降低估计的协方差。
3.5.3倒频谱分析：
图3-8实倒谱
图3-9复倒谱
正弦信号，其第一个功率谱变换为一脉冲，经滤波后进入第二次功率谱变换，其输出为幅度很低的三角波输出，因而检测不到其存在。
3.5.4细化谱分析：
图3-10原始信号FFT
图3-11ZOOM-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）;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（'原始信号'）;
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）
subplot（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;%采样点数为1024
Pxx_short=abs（fft（x,256））.^2/256;%采样点数为256
figure（3）
subplot（2,1,1）
plot（（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+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（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-FFT
fs=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;
loadzoomfftdata;
fi=6;%最小细化截止频率
np=10;%放大倍数
nfft=512;%fft长度
nt=length（x）;
fa=fi+0.5*fs/np;%最大细化截止频率
nf=2^nextpow2（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）*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（nn））,':
',f2,abs（y2））;
%存储ZOOM-FFT结果
saveafterzoomdataf2y2
学习心得
通过这将近一周半的数字信号处理的课程设计，我先在图书馆里查找了相关的书籍，如MATLAB类的编程书籍，各类数据处理类的书籍以及机械振动类的书籍等，即丰富了自己的知识范围，又对与自己所学的知识有了更深的了解和认识，同时也对它的应用有了一个大体的认识。
这样将会更加激励我好好学习相关的知识，不断的将所学的知识用于实践。
于实践中牢牢的掌握它。
在设计的过程中，我也认识到了自己所学知识的不足。
这也让我再次认识到知识是无尽的，只有不断的充实自己、完善自己的知识理论体系，才能够更好的胜任自己以后的工作。
设计过程中知识的不足也让我更加坚定了终身学习的决心。
在设计的过程中，我也得到了我们设计小组的成员和很多同学的帮组。
这也加强了我与其他同学合作的能力。
查找资料的过程中我也增强自己学习的能力，这些都将让我在以后的学习、生活和工作中受益匪浅。
总之，对于这样的课程设计活动，我收获了很多东西，也将使我在以后的学习、工作中更加轻松和积极。
这也正是参加这次活动的目的和意义。
参考文献
1.林洪彬，谢平，王娜.信号处理原理及应用.机械工业出版社，2009年．
2.李力．机械信号处理及其应用．华中科技大学出版社，2007年．
3.罗军辉，白义臣．MATLAB7.0在数字信号处理中的应用．机械工业出版社，2005年．
4.周浩敏，王睿.测试信号处理技术.北京航空航天大学出版社，2005年．
5.程佩青.数字信号处理教程,清华大学，2001年.
6.邓立新.数字信号处理学习辅导及习题详解,电子工业，2003年.
7.丁玉美，数字信号处理（2ed）学习指导,西安电子科技大学,2001年.
8.胡广书，数字信号处理-理论算法与实现,清华大学，2003年.
升级会员 