机械动态信号分析与处理提交结题参考

东北大学
研究生考试试卷
考试科目：机械动态信号分析与处理课程编号：09581066
阅卷人：
考试时间：2013.1.9
姓名：叶利丽
学号：1370285
东北大学研究生院
Matlab程序报告
一、设计题目
用MATLAB编写程序对机械振动信号进行分析
二、设计目的
去除信号中混杂的白噪声并对其在时、频域上进行分析。

三、设计的理论基础
(1) MATLAB的简单介绍
MATLAB是美国Mathworks公司开发的新一代科学计算软件：MATLAB是英文MATtrix LABoratory(矩阵实验室)的缩写；MATLAB是一个专门为科学计算而设计的可视化计算器。

利用这个计算器中的简单命令，能快速完成其他高级语言只有通过复杂此案出才能实现的数值计算和图形显示。

MATLAB是一种既可交互使用又能解释执行的计算机编程语言。

所谓交互使用，是指用户输入一条语句后立即就能得到该语句的计算结果，而无需像C语言那样首先编写源程序，然后对之进行编译，连接，才能最终形成可执行文件。

MATLAB语言可以用直观的数学表达式来描述问题，从而避开繁琐的底层编程，因此可大大提高工作效率。

MATLAB是解决工程技术问题的技术平台。

利用它能够轻松完成复杂的数值计算，数据分析，符号计算和数据可视化等任务。

MATLAB软件由主包和各类工具箱构成。

其中，主包基本是一个用C/C++等语言编写成的函数库。

该函数库提供矩阵（或数组）的各种算法以及建立在此基础上的各种应用函数和一些相关的用户有好操作界面。

而工具箱从深度和广度上大大扩展了MATLAB主包的功能和应用领域。

随着自身的不断完善和发展，MATLAB功能越来越强大，应用也越来越广泛。

(2) 信号测试技术与分析
随着机械工业不断向自动化、高精度、智能化等方向的发展，在机械设备运行及生产过程中进行参量测试、分析与诊断等处理过程已成为必要环节，许多信号处理方法如时域统计分析、相关分析、相干分析、频谱分析等已经被广泛被应用与机械工程测试领域。

测试信号通常指的是被测对象的运动或状态信息。

测试信号可以用数学表达
式描述，也可以用图形、图表等进行描述。

在工程测试中，有的信号可以用数学公式精确描述，而大量的测试信号却只能用数学公式来近似描述。

随着数字信号处理理论和快速算法的不断完善，特别是数字信号处理芯片的高速发展，测试信号分析与处理技术已经在许多领域越来越普及。

本设计为机械测试信号的时域、频域分析以及降噪，其中时域分析包括对信号最大值、最小值、中值、方差、相关的分析，频域分析包括FFT分析、频谱分析、功率谱分析、降噪主要包括滤波、小波分析以及小波重构等。

时域统计分析是指对动态信号的各种时域参数、指标的估计或运算，通过选择和考察合适的信号动态分析指标，可以对不同类型的故障做出准确的判断。

相关分析方法是对机械信号进行时域分析方法的常用之一，也是故障诊断的重要手段。

无论是分析两个随机变量之间的关系，还是分析两个信号或一个信号在一定时移前后之间的关系，都需要应用相关分析。

所谓相关，就是指变量之间的线性联系或相互依赖关系。

频域分析中FFT是离散傅立叶变换DFT的快速算法，它使DFT在数字通信、语音信号处理、图像处理、系统仿真、雷达理论、光学、医学、地震监测等各个领域得到了广泛应用。

然而，各种应用一般都以卷积和相关运算的具体处理为依据，或者以DFT作为连续傅立叶变换的基础。

FFT主要包括两类：按时间抽取法和按频率抽取法。

倒频谱实际上是频域信号取对数的傅立叶变换再处理，或称为“频域信号的傅立叶变换再变换”。

对功率谱密度函数取对数的目的是使再变换以后信号的能量更加集中。

倒频谱可以分析复杂频谱上的周期成分，分离和提取在密集泛频信号中的成分。

对于具有同族谐频和异族谐频等复杂信号的分析，效果很好。

功率谱估计是数字信号处理的主要内容之一，主要研究信号在频域中的各种特征，目的是根据有限数据在频域内提取被淹没在噪声中的有用信号。

滤波是根据信号的不同特性的消除或减弱干扰噪声，提取有用信号的过程，对数字信号的滤波主要用FIR数字滤波器。

FIR滤波器总是稳定的，能够用硬件实现，FIR滤波器具有线性相位。

小波分解是对信号采样后，可得到在一个大的有限频带中的一个信号，对这个信号进行小波多尺度分解，其实质就是把采到的信号分成两个信号，即高频部
分和低频部分，而低频部分通常包含了信号的主要信息，高频部分则与噪音及扰动联系在一起。

根据分析的需要，可以继续对所得到的低频部分进行分解，如此又得到了更低频部分的信号和频率相对较高部分的信号。

四、基础函数
构造一个振动信号，有三种频率和白噪声组成的信号，再利用MATLAB中的函数mean()、min()、max()、std()、xcorr()对离散序列中的平均值、最大值、最小值、标准差、自相关等时域性能进行分析，通过调用函数fft()、power()；对该振动信号进行频域内的性能分析。

用filter()、wavedec()、detcoef()、appcoef()、wavere()进行滤波、小波分解与小波重构。

五、实现过程
(1) 原始信号
fs=1000;
t=0:0.001:1;
x=3*sin(2*pi*10*t)+2*sin(2*pi*20*t)+5*sin(2*pi*30*t)+randn(size(t));
figure(1);plot(x);
title('振动信号的波形');
xlabel('时间/秒');
ylabel('幅度');
grid; hold on;
图(1)为原始信号的波形，他由三个正弦波和一个白噪声信号组成。

信号中有周期成分并且有小毛刺。

020*******
80010001200
-15-10
-5
5
10
15
振动信号的波形
时间/秒幅度
图(1) 原始信号的波形
（2）时域分析
时域分析结果
序列的平均值为 -0.0319
序列的最小值为 -10.3266
序列的最大值为 11.7024
序列的标准差为 4.4622
图(2)为信号自相关函数的波形，从图中可以看出信号中当 很大时信号逐渐衰减，说明信号中含有非周期成分，此成分为白噪声引起的。

图(3)为自相关函数的功率谱，从中可以看出所对应的周期信号的频率。

-1-0.8-0.6-0.4-0.20
0.20.40.60.81
-1.5-1
-0.5
0.5
1
1.5
2
4自相关函数图
时间/s R x (t )
图(2)信号的自相关函数
050100150200250300
350400450500
01000
2000
3000
4000
5000
6000
7000
自相关频谱
频率/Hz
图(3) 自相关函数的频谱
(3)频域分析
图(4)为信号的频谱，从信号的频谱中可以看出信号在10HZ 、20HZ 、30HZ 时有峰值，且此时对应的幅值分别接近3、2、5。

信号中也存在着微小的振动量为白噪声。

050100150200250300
350400450500
00.5
1
1.5
2
2.5
3
3.5
4
4.5
5
振动信号频谱
频率/HZ 幅值图(4) 信号的频谱
050100150200250300
350400450500
05
10
15
20
25
振动信号功率谱
频率/HZ 幅值.2
图(5)振动信号的功率谱
（5）滤波
采用FIR 滤波器进行信号的50阶带通滤波，滤波的带宽为[10,30]。

图(6)为滤波后的信号。

信号滤波后还是很干净的。

00.10.20.30.40.5
0.60.70.80.91
-8-6
-4
-2
2
4
6
8
10
滤波后的信号
时间/s 幅度
图(6)滤波后的信号
下面用小波分析的方法进行的降噪处理，图(7)为小波分解后的函数。

图(8)对小波函数的的重构。

020406080100120140
-50
50
c a 3小波分解
020406080100120140
-10
10c d 3050100150200250300
-5
5c d 20100200300400500600-5
5c d 1020040060080010001200-200
20
x
图(7) 小波分解后的信号
020040060080010001200
-10-8
-6
-4
-2
2
4
6
8
10
去噪之后的波形
图(8) 小波重构信号
附件：程序代码
%%三种信号和一个白噪声的信号分析并去噪%%
close all;
clear all;
%原始信号%
fs=1000;
t=0:0.001:1;
x=3*sin(2*pi*10*t)+2*sin(2*pi*20*t)+5*sin(2*pi*30*t)+randn(size(t)); figure(1);plot(x);
title('振动信号的波形');
xlabel('时间/秒');
ylabel('幅度');
grid; hold on;
%时域分析%
m= mean(x); disp (m);
mi= min(x); disp (mi);
mx = max(x); disp (mx);
st = std(x); disp (st);
%自相关分析%
[w,lags]=xcorr(x);
figure(7);
plot(lags/fs,w);
title('自相关函数图');
xlabel('时间/s');
ylabel('Rx(t)');
grid; hold on;
%自相关功率谱%
m=length(w)-1;
z=fft(w,m);
mag_z=abs(z);
ff=(0:m-1)*fs/m;
figure(8);
plot(ff(1:m/2),mag_z(1:m/2)*2/m); title('自相关频谱'),xlabel('频率/Hz'); %频域分析%
l=length(x);
r=fft(x);
r=1/500*abs(r);
r=r(1:(end-1)/2);
f=linspace(0,0.5,(l-1)/2)*fs;
figure(2);plot(f,r);
title('振动信号频谱');
xlabel('频率/HZ');
ylabel('幅值');
grid; hold on;
%功率谱%
figure(3);
power=r.^2;
plot(f,power);
title('振动信号功率谱');
xlabel('频率/HZ');
ylabel('幅值.^2');
grid; hold on;
%滤波%
wn=[10 30]*2/fs;
b=fir1(40,wn);
y=filter(b,1,x);
figure(4);
plot(t,y);
title('滤波后的信号');
xlabel('时间/s');
ylabel('幅度');
grid; hold on;
%小波分析%
[c,l]=wavedec(x,3,'db4');
[cd1,cd2,cd3]=detcoef(c,l,[1,2,3]);
ca3=appcoef(c,l,'db4',3);
figure(5);
subplot(511);plot(1:l(1),ca3);ylabel('ca3');title('小波分解'); subplot(512);plot(1:l(2),cd3);ylabel('cd3');
subplot(513);plot(1:l(3),cd2);ylabel('cd2');
subplot(514);plot(1:l(4),cd1);ylabel('cd1');
subplot(515);plot(x);ylabel('x');
grid; hold on;
%小波重构%
figure(6);
cdd3=zeros(1,length(cd3));
cdd2=zeros(1,length(cd2));
cdd1=zeros(1,length(cd1));
c1=[ca3 cdd3 cdd2 cdd1];
s1=waverec(c1,l,'db4');
plot(s1);
title('去噪之后的波形');
grid; hold on;。