南昌航空大学实验报告课程名称：数字信号处理实验名称：滚动轴承的振动信号特征分析实验时间： 2013年5月14日班级： 100421学号： **********：***成绩：滚动轴承的振动信号特征分析一、实验目的利用《数字信号处理》课程中学习的序列运算、周期信号知识、DFT 知识，对给定的正常轴承数据、内圈故障轴承数据、外圈故障轴承数据、滚珠故障轴承数据进行时域特征或频域特征提取和分析，找出能区分四种状态（滚动轴承的外圈故障、内圈故障、滚珠故障和正常状态）的特征。

二、实验原理振动机理分析：机械在运动时，由于旋转件的不平衡、负载的不均匀、结构刚度的各向异性、间隙、润滑不良、支撑松动等因素，总是伴随着各种振动。

振动的幅值、频率和相位是振动的三个基本参数，称为振动三要素。

幅值：幅值是振动强度的标志，它可以用峰值、有效值、平均值等方法来表示。

频率：不同的频率成分反映系统内不同的振源。

通过频谱分析可以确定主要频率成分及其幅值大小，从而寻找振源，采取相应的措施。

相位：振动信号的相位信息十分重要，如利用相位关系确定共振点、测量振型、旋转件动平衡、有源振动控制、降噪等。

对于复杂振动的波形分析，各谐波的相位关系是不可缺少的。

在振动测量时，应合理选择测量参数，如振动位移是研究强度和变形的重要依据；振动加速度与作用力或载荷成正比，是研究动力强度和疲劳的重要依据；振动速度决定了噪声的高低，人对机械振动的敏感程度在很大频率范围内是由速度决定的。

速度又与能量和功率有关，并决定动量的大小。

提取振动信号的幅域、时域、频域、时频域特征，根据特征进行故障有无、故障类型和故障程度三个层次的判断。

三、 实验内容Step1、使用importdata （）函数导入振动数据。

Step2、把大量数据分割成周期为单元的数据，分割方法为：设振动信号为{x k }(k =1,2,3,…,n )采样频率为f s ，传动轴的转动速率为V r 。

采样间隔为：1st f ∆= （1）旋转频率为：60rr V f =（2）传动轴的转动周期为：1r T f = （3）由式（1）和（3）可推出振动信号一个周期内采样点数N ：11s r r sf f TN t f f ===∆ （4） 由式（2）可得到传动轴的转动基频f r =29.95Hz ，再由式（3）可得到一个周期内采样点数N=400.67，取N =400。

Step3、提取振动信号的特征，分析方法包括：1、时域统计分析指标（波形指标(Shape Factor)、峰值指标(Crest Factor)、脉冲指标(Impulse Factor)、裕度指标(Clearance Factor)、峭度指标(KurtosisValue) ）等，相关计算公式如下：（1）波形指标：Pf X WK X =(5)其中，PX 为峰值，X 为均值。

p X 计算公式如下：),,2,1}(max{N i x X i p ==（2）波峰指标：pf rms X CK X =(6)其中，rmsX 为均方根值。

均方根值计算如下：rms X （3）峭度指标：4rmsK X β=(8)其中，β为峭度，计算式如下：411Ni i x N β==∑ (9)（4）脉冲指标：PX I X =(10)（5）偏度指标：3rmsP X α=(11)其中，α为偏度，计算式如下：311Nii xNα==∑ (12)（6）裕度指标：rp X X L =(13)其中，r X 为方根幅值，计算公式如下：211⎥⎦⎤⎢⎣⎡=∑=N i i r x N X2、频域分析，包括离散频谱序列统计分析指标，阶次分析等。

滚动轴承的理论特征频率计算公式如下： （1）内圈故障频率： (1cos())2i r N df f Dα=+ =161.73 Hz （2）外圈故障频率： (1cos())2o r N df f Dα=-=107.82 Hz （3）滚珠故障频率： 22(1()cos ())2b r D df f d Dα=-=141.17 Hz 其中，N 为滚珠数，d 为滚珠直径，α 为接触角,D 为轴承节径。

Step4（提高要求）、提取未知信号的特征，与已知信号的特征进行比对分类，诊断信号是否存在故障以及故障类型。

未知故障类型信号：1.mat四、实验的程序代码（1）、待测信号1function [y1,y2,y3,y4]=KuttosisFactor(); s=importdata('E:\原始信号\内圈故障\105.mat'); s_value=s.X105_DE_time'; s_cutvalue=s_value(1:24000); tree=wpdec(s_cutvalue,3,'sym6'); det=wpcoef(tree,2);sigma=median(abs(det))/0.675;alpha=2;thr=wpbmpen(tree,sigma,alpha);keepapp=1;xd=wpdencmp(tree,'s','nobest',thr,keepapp);T=400;Vaule1=zeros(1,50);for n=1:1:50B=sum(xd((1+(n-1)*T):(n*T)).^4)/T;Rms=sqrt(sum(xd((1+(n-1)*T):(n*T)).^2)/T); Vaule1(n)=B/(Rms^4);end;s=importdata('E:\原始信号\正常\97.mat');s_value1=s.X097_DE_time';s_cutvalue=s_value1(1:24000);tree=wpdec(s_cutvalue,3,'sym6');det=wpcoef(tree,2);sigma=median(abs(det))/0.675;alpha=2;thr=wpbmpen(tree,sigma,alpha);keepapp=1;xd=wpdencmp(tree,'s','nobest',thr,keepapp); Vaule2=zeros(1,50);for n=1:1:50B=sum(xd((1+(n-1)*T):(n*T)).^4)/T;Rms=sqrt(sum(xd((1+(n-1)*T):(n*T)).^2)/T); Vaule2(n)=B/(Rms^4);end;s=importdata('E:\原始信号\外圈故障\130.mat'); s_value1=s.X130_DE_time';s_cutvalue=s_value1(1:24000);tree=wpdec(s_cutvalue,3,'sym6');det=wpcoef(tree,2);sigma=median(abs(det))/0.675;alpha=2;thr=wpbmpen(tree,sigma,alpha);keepapp=1;xd=wpdencmp(tree,'s','nobest',thr,keepapp); Vaule3=zeros(1,50);for n=1:1:50B=sum(xd((1+(n-1)*T):(n*T)).^4)/T;Rms=sqrt(sum(xd((1+(n-1)*T):(n*T)).^2)/T); Vaule3(n)=B/(Rms^4);end;s=importdata('E:\原始信号\滚珠故障\118.mat'); s_value1=s.X118_DE_time';s_cutvalue=s_value1(1:24000);tree=wpdec(s_cutvalue,3,'sym6');det=wpcoef(tree,2);sigma=median(abs(det))/0.675;alpha=2;thr=wpbmpen(tree,sigma,alpha);keepapp=1;xd=wpdencmp(tree,'s','nobest',thr,keepapp); Vaule4=zeros(1,50);for n=1:1:50B=sum(xd((1+(n-1)*T):(n*T)).^4)/T;Rms=sqrt(sum(xd((1+(n-1)*T):(n*T)).^2)/T); Vaule4(n)=B/(Rms^4);end;s=importdata('E:\待诊断信号\1.mat');s_value1=s;s_cutvalue=s_value1(1:24000);tree=wpdec(s_cutvalue,3,'sym6');det=wpcoef(tree,2);sigma=median(abs(det))/0.675;alpha=2;thr=wpbmpen(tree,sigma,alpha);keepapp=1;xd=wpdencmp(tree,'s','nobest',thr,keepapp);T=400;Vaule5=zeros(1,50);for n=1:1:50B=sum(xd((1+(n-1)*T):(n*T)).^4)/T;Rms=sqrt(sum(xd((1+(n-1)*T):(n*T)).^2)/T);Vaule5(n)=B/(Rms^4);end;y1(1)=(Vaule1(1)+Vaule1(2)+Vaule1(3))/3;y2(1)=(Vaule2(1)+Vaule2(2)+Vaule2(3))/3;y3(1)=(Vaule3(1)+Vaule3(2)+Vaule3(3))/3;y4(1)=(Vaule4(1)+Vaule4(2)+Vaule4(3))/3;y5(1)=(Vaule5(1)+Vaule5(2)+Vaule5(3))/3;y1(50)=(Vaule1(50)+Vaule1(49)+Vaule1(48))/3;y2(50)=(Vaule2(50)+Vaule2(49)+Vaule2(48))/3;y3(50)=(Vaule3(50)+Vaule3(49)+Vaule3(48))/3;y4(50)=(Vaule4(50)+Vaule4(49)+Vaule4(48))/3;y5(50)=(Vaule5(50)+Vaule5(49)+Vaule5(48))/3;for n=2:1:49y1(n)=(Vaule1(n-1)+Vaule1(n)+Vaule1(n+1))/3;y2(n)=(Vaule2(n-1)+Vaule2(n)+Vaule2(n+1))/3;y3(n)=(Vaule3(n-1)+Vaule3(n)+Vaule3(n+1))/3;y4(n)=(Vaule4(n-1)+Vaule4(n)+Vaule4(n+1))/3;y5(n)=(Vaule5(n-1)+Vaule5(n)+Vaule5(n+1))/3;end;%绘图k=linspace(1,50,50);plot(k,y1(k),'--.k',k,y2(k),'-r',k,y3(k),'-pb',k,y4(k),'--*g',k,y5(k),'m:d');legend('内圈故障','正常''外圈故障',’故障’,’诊断信号’);小结：通过现象观察可以看出，以上得到的图形都是外圈故障。