傅里叶分析滚动轴承的故障诊断

作业名称：傅里叶分析滚动轴承的故障诊断院系：机械工程系学号：姓名：指导教师：20XX年XX月XXXXXXXXX校区傅里叶分析滚动轴承的故障诊断摘要：简要介绍了快速傅里叶变换（FFT）在滚动轴承故障分析中的应用，滚动轴承在机械设备中使用非常广泛，其工作状态直接影响整个设备的运行品质。对滚动轴承进行状态监测与故障诊断，能够避免重大事故的发生，获得较大的经济和社会效益。通过快速傅里叶变换（FFT）对滚动轴承运行时的实时数据信号进行分析，可以实现对滚动轴承的状态监测和故障诊断。同时，采用对正常轴承和故障轴承信号对比分析、各种故障轴承之间信号的对比分析，加深了快速傅里叶变换（FFT）对轴承实时信号分析的运用和理解，能够更好的对轴承进行状态监测和故障分析。关键词：快速傅里叶变换（FFT）；滚动轴承；故障诊断；状态监测Abstract:ThispaperdescribesafastFouriertransform(FFT)intherollingbearingfailureanalysisapplications,bearinginmachineryandequipmentiswidelyused,anditsworkingstatusdirectlyaffectsthequalityoftheoperationoftheentiredevice.Rollingelementbearingconditionmonitoringandfaultdiagnosis,abletoavoidmajoraccidentsandachievegreatereconomicandsocialbenefits.ThroughFastFourierTransform(FFT)forreal-timedatabearingsignalruntimeanalysiscanbeachievedontherollingbearingconditionmonitoringandfaultdiagnosis.Meanwhile,theuseofnormalbearingsandbearingfaultsignalcomparativeanalysisofvariousfaultsignalscomparativeanalysisbetweenthebearingsanddeepenedthefastFouriertransform(FFT)ofthebearingusingreal-timesignalanalysisandunderstandingofthebearingcanbebetterconditionmonitoringandfaultanalysis.Keywords:fastFouriertransform(FFT);Rolling;faultdiagnosis;conditionmonitoring一、概述通过对快速傅里叶变换（FFT）的原理的理解和学习，利用MATLAB软件编程应用快速傅里叶变换（FFT）的方法，对滚动轴承的1组正常数据和2组故障数据（故障类型不同）进行信号分析和处理，并对正常轴承和故障轴承信号对比分析、各种故障轴承之间信号的对比分析，并得出结论，实现对滚动轴承的状态监测和故障分析。二、信号处理方法及原理快速傅里叶变换，是计算离散傅里叶变换（DFT）的一种快速算法,简称FFT。当用数字计算机计算信号序列x(n)的离散傅里叶变换时，它的正变换(1)反变换(IDFT)是(2)式中、x(n)和X(k)可以是实数或复数。由上式可见,要计算一个抽样序列就需要做N次复数乘法运算及N-1次复数加法运算。计算离散傅里叶变换的快速方法，有按时间抽取的FFT算法和按频率抽取的FFT算法。前者是将时域信号序列按偶奇分排,后者是将频域信号序列按偶奇分排。它们都借助于的两个特点:一是的周期性；另一是的对称性，这里符号*代表其共轭。这样,便可以把离散傅里叶变换的计算分成若干步进行，计算效率大为提高。时间抽取算法令信号序列的长度为N＝2M,其中M是正整数，可以将时域信号序列x(n)分解成两部分，一是偶数部分x（2n），另一是奇数部分x（2n+1），其中。于是信号序列x(n)的离散傅里叶变换可以用两个N/2抽样点的离散傅里叶变换来表示和计算。考虑到和离散傅里叶变换的周期性,式(1)可以写成(3)其中(4a)(4b)由此可见，式(4)是两个只含有N/2个点的离散傅里叶变换，G(k)仅包括原信号序列中的偶数点序列，H(k)则仅包括它的奇数点序列。虽然k=0，1，2，…，N-1，但是G(k)和H(k)的周期都是N/2,它们的数值以N/2周期重复。因为于是由式(3)和式(4)得到(5a)(5b)因此,一个抽样点数为N的信号序列x(n)的离散傅里叶变换，可以由两个N/2抽样点序列的离散傅里叶变换求出。依此类推，这种按时间抽取算法是将输入信号序列分成越来越小的子序列进行离散傅里叶变换计算，最后合成为N点的离散傅里叶变换。通常用图1中蝶形算法的信号流图来表示式(5)的离散傅里叶变换运算。例如，N＝8＝23的抽样点的信号序列x(n)的离散傅里叶变换,可用如图2所示的FET算法的信号流图来计算。由图可知：①N=2M点的离散傅里叶变换的计算全由蝶形运算组成，需要M级运算,每级包括N/2个蝶形运算，总共有kK个蝶形运算。所以，总的计算量为次复数乘法运算和Nlog2N次复数加法运算。②FFT算法按级迭代进行，计算公式可以写成(6)N抽样点的输入信号具有N个原始数据x0(n),经第一级运算后，得出新的N个数据x1(n)，再经过第二级迭代运算，又得到另外N个数据x2(n)，依此类推，直至最后的结果x(k)＝xM(k)＝X(k)在逐级迭代计算中，每个蝶形运算的输出数据存放在原来存贮输入数据的单元中，实行所谓“即位计算”，这样可以节省大量存放中间数据的寄存器。③蝶形运算中加权系数随迭代级数成倍增加。由图2可以看出系数的变化规律。对于N=8,M=3情况,需进行三级迭代运算。在第一级迭代中，只用到一种加权系数;蝶形运算的跨度间隔等于1。在第二级迭代中，用到两种加权系数即、；蝶形运算的跨度间隔等于2。在第三级迭代中，用到4种不同的加权系数即、、、；蝶形运算的跨度间隔等于4。可见,每级迭代的不同加权系数的数目比前一级迭代增加一倍；跨度间隔也增大一倍。④输入数据序列x(n)需重新排列为x(0)、x(4)、x(2)、x(6)、x(1)、x(5)、x(3)、x(7),这是按照二进制数的码位倒置所得到的反序数，例如N=8中数“1”的二进制数为“001”，将其码位倒转变为“100”，即为十进制数“4”。频率抽取算法按频率抽取的FFT算法是将频域信号序列X(k)分解为奇偶两部分，但算法仍是由时域信号序列开始逐级运算，同样是把N点分成N/2点计算FFT，可以把直接计算离散傅里叶变换所需的N2次乘法缩减到次。在N＝2的情况下,把N点输入序列x(n)分成前后两半(7)时间序列x1(n)±x2(n)的长度为N/2,于是N点的离散傅里叶变换可以写成(8a)(8b)频率信号序列X(2l)是时间信号序列x1(n)+x2(n)的N/2点离散傅里叶变换,频率信号序列X(2l+1)是时间信号序列【x1(n)-x2(n)】的N/2点离散傅里叶变换,因此,N点离散傅里叶变换的计算，通过两次加（减）法和一次乘法，从原来序列获得两个子序列，所以，频率抽取算法也具有蝶形运算形式。以2为基数的FFT基本蝶形运算公式为(9)其计算量完全和时间抽取算法一样，即只需次乘法运算和Nlog2N次加(减)法运算。图3表示N=8=23点的离散傅里叶变换的信号流图。由图可见,它以三级迭代进行即位计算,输入数据是按自然次序存放，使用的系数也是按自然次序，而最后结果则以二进制反序存放。实际上，频率抽取算法与时间抽取算法的信号流图之间存在着转置关系，如将流图适当变形，可以得出多种几何形状。除了基2的FFT算法之外，还有基4、基8等高基数的FFT算法以及任意数为基数的FFT算法。三、故障诊断的结果选取正常轴承数据normal2.mat，内圈故障数据inner-race2.mat，外圈故障数据outer-race2.mat，进行数据信号分析，得出信号时域图和信号频谱图。分别如图4、图5和图6所示。图4.normal2.mat处理结果图5.inner-race2.mat处理结果图6.outer-race2.mat处理结果从正常轴承的频谱图(图4)可以看出，在频率为0~2000Hz和10000~12000Hz的频段有较高阶谐波，且呈对称状态，幅值较大，最大幅值在1000Hz和11000Hz左右。在2000~10000Hz的频段中，幅值很小。从内圈故障的频谱图(图5)可以看出，在频率为0~4000Hz和8000~12000Hz的频段有较高阶谐波，且呈对称状态。在4000~8000Hz的频段中，波形幅值较小。从外圈故障的频谱图(图6)可以看出，在频率为0~5000Hz和7000~12000Hz的频段有较高阶谐波，且呈对称状态，最大幅值在1000Hz和11000Hz左右。在5000~7000Hz的频段中，波形振幅较小。四、结论通过此次对滚动轴承的故障检测和分析，使我获益良多。但是由于各种特征频率都是从理论上推导出来的，而实际上，由于轴承的各几何尺寸会有误差，加上轴承安装后的变形、FFT计算误差等因素，使得实际的频率与计算所得的频率会有些出入。所以在频谱图上寻找各特征频率时，须在计算的频率值上找其近似值来作诊断。通过此次学习，加深了我对MATLAB的熟悉，使我更加熟练掌握了MATLAB，同时更加理解和掌握了FFT的原理和方法。附：MATLAB程序（1）正常轴承程序x=X098_DE_time;%信号数组subplot(2,1,1);plot(x);%时域波形xlabel('时间序列');ylabel('幅值');title('信号时域图');fs=12000;%采样频率N=length(x);n=0:N-1;y=fft(x,N);%进行fft变换m=abs(y(1:N))*2/N;%求信号的真实幅值f=n*fs/N;%进行对应的频率转换subplot(2,1,2)stem(f(1:N),m(1:N));%绘出频谱图xlabel('频率/Hz');ylabel('幅值');title('信号频谱图');gridon;（2）内圈故障轴承程序x=X274_DE_time;%信号数组subplot(2,1,1);plot(x);%时域波形xlabel('时间序列');ylabel('幅值');title('信号时域图');fs=12000;%采样频率N=length(x);n=0:N-1;y=fft(x,N);%进行fft变换m=abs(y(1:N))*2/N;%求信号的真实幅值f=n*fs/N;%进行对应的频率转换subplot(2,1,2)stem(f(1:N),m(1:N));%绘出频谱图xlabel('频率/Hz');ylabel('幅值');title('信号频谱图');gridon;（3）外圈故障轴承程序x=X313_DE_time;%信号数组subplot(2,1,1);plot(x);%时域波形xlabel('时间序列');ylabel('幅值');title('信号时域图');fs=12000;%采样频率N=length(x);n=0:N-1;y=fft(x,N);%进行fft变换m=abs(y(1:N))*2/N;%求信号的真实幅值f=n*fs/N;%进行对应的频率转换subplot(2,1,2)stem(f(1:N),m(1:N));%绘出频谱图xlabel('频率/Hz');ylabel('幅值');title('信号频谱图');gridon;