數(shù)學(xué)建模：齒輪箱故障診斷

表單

gearbox00

為齒輪箱正常工況下采集到的振動信號;表單

gearbox10

為故障狀態(tài)

1

下采集到的振動信號;表單

gearbox20

為故障狀態(tài)

2

下采集到的故障信號;表單

gearbox30

為故障狀態(tài)

3

下采集到的故障信號;表單

gearbox40

為故障狀態(tài)

4

下采集到的振動信號。1、對齒輪箱各個狀態(tài)下的振動數(shù)據(jù)進行分析，研究正常和不同故障狀態(tài)下振動數(shù)據(jù)的變化規(guī)律及差異，并給出刻畫這些差異的關(guān)鍵特征。這題每個狀態(tài)有四個指標(biāo)，即對每個狀態(tài)的各指標(biāo)首先做圖觀察，看看每個狀態(tài)下的數(shù)據(jù)變化趨勢，明天我簡單做做給出部分圖和代碼，不過用spss或者Excel也是一樣可以的。至于差異的關(guān)鍵特征，對數(shù)據(jù)進行特征分析即可。代碼也簡單，主要就是看寫論文吧。2、建立齒輪箱的故障檢測模型，對其是否處于故障狀態(tài)進行檢測，并對模型的性能進行評價。

根據(jù)第一問提取的數(shù)據(jù)特征通過計算故障強度系數(shù)再進行后續(xù)計算，見下文《普通公路隧道機電設(shè)施故障檢測方法研究》用小波分析法檢測故障，見下文《基于向量自回歸模型和小波分析法的列車充電機電流傳感器故障檢測方法》3、建立齒輪箱的故障診斷模型，對其處于何種故障狀態(tài)進行判斷，并對模型的性能進行評價。根據(jù)前一問的結(jié)果先檢測哪些是故障，后用一種基于樹形卷積神經(jīng)網(wǎng)絡(luò)模型診斷是哪種故障。見下文《鐵路信號聯(lián)鎖故障診斷模型構(gòu)建及仿真》4、結(jié)合所建立的故障檢測和診斷模型對附件

2

中另行采集的

12

組測試數(shù)據(jù)進行檢測和診斷分析，將分析結(jié)果填寫到下表中(注:測試數(shù)據(jù)中可能存在除以上

4

種故障之外的故障狀態(tài)，若存在，則將對應(yīng)的診斷結(jié)果標(biāo)記為:其它故障)，并將此表格放到論文的正文中第四問直接運用前兩問的模型即可。明天繼續(xù)更新，今天時間有點倉促，只能大概提供這一點點。明天我會進行修改和補充，大家有問題評論或私信即可，明天我的課還是比較多的。所以應(yīng)該周末才會給出詳細的思路和部分代碼。一樣的，對于數(shù)據(jù)題我們首先要做的是數(shù)據(jù)預(yù)處理，接下來再進行下面的操作。關(guān)于問題一：

可以做每組數(shù)據(jù)的散點圖，我這個就給出一張例圖。Excel和SPSS都可以做的，SPSS做的估計會更好看一點。Excel好好調(diào)參也可以很好看。做這種圖的，第一自己觀察看看有沒有什么顯著差異，第二寫論文。這個也算是震動信號數(shù)據(jù)的一種分析，論文寫寫就行。后面的計算有比較多的方法1、通過計算方差，篩選特征。計算方法參考下文\o"??????數(shù)據(jù)篩選特征方法-方差法_gao_vip的博客-CSDN博客_方差選擇法特征篩選"??????數(shù)據(jù)篩選特征方法-方差法_gao_vip的博客-CSDN博客_方差選擇法特征篩選2、計算特征重要性，取重要性高的那類作為特征。3、SVD（奇異值分解），這個方法比較好的。大家可以去看看下面這個論文，這一問能搞定?！禨VD曲率譜降噪和快速譜峭度的滾動軸承微弱故障特征提取》\o"常用數(shù)據(jù)特征提取，時域特征、頻域特征、小波特征提取匯總；特征提??；有效matlab代碼_入間同學(xué)的博客-CSDN博客_小波熵特征提取"常用數(shù)據(jù)特征提取，時域特征、頻域特征、小波特征提取匯總；特征提??；有效matlab代碼_入間同學(xué)的博客-CSDN博客_小波熵特征提取4、小波分析特征提取%裝入變換放大器輸入輸出數(shù)據(jù)%bf_150ms.dat為正常系統(tǒng)輸出信號%bf_160ms.dat為故障系統(tǒng)輸出信號loadbf_150ms.dat;loadbf_160ms.dat;s1＝bf_150ms(1:1000);%s1為正常信號s2＝bf_160ms(1:1000);%s2為故障信號%畫出正常信號與故障信號的原始波形tittle(“原始信號’);Ylabel('s1');subplot(922);plot(s2);title('故障信號');Ylabel('s2');%============================================%用dbl小波包對正常信號s1進行三層分解[t,d]＝wpdec(sl,3,'db','shannon');%plontree(t)%畫小波包樹結(jié)構(gòu)的圖形%下面對正常信號第三層各系數(shù)進行重構(gòu)%s130是指信號sl的[3,0]結(jié)點的重構(gòu)系數(shù);其他依次類推sl30＝wprcoef(t,d,[3,0]);s13l＝wprcoef(t,d,[3,1]);s132＝wprcoef(t,d,[3,2]);sl33＝wprcoef(t,d,[3,3]);sl34＝wprcoef(t,d,[3,4]);s135＝wprcoef(t,d,[3,5]);s136＝wprcoef(t,d,[3,6]);s137＝wprcoef(t,d,[3,7]);%畫出至構(gòu)系數(shù)的波形subplot(9,2,3);plot(s130);Ylabel('S130');subpolt(9,2,5);plot(s131);Ylabel('S13l');subplot(9,2,7);plot(s132);Ylabel('S132');subplot(9,2,9);plot(s133);Ylabel('S133');subplot(9,2,11);plot(s134);Ylabel('S134');subplot(9,2,13);plot(s135);Ylabel('S135');subplot(9,2,15);plot(s136);Ylabel('S136');subplot(9,2,17);plot(s137);Ylabel('S137');%--------------------------------------%計算正常信號各重構(gòu)系數(shù)的方差%s10是指s130的方差，其他依此類推s10＝norm(sl30);sll＝norm(s131);s12＝norm(sl32);s13＝norm(sl33);sl4＝norm(s134);s15＝norm(s135);s16＝norm(sl36);s17＝norm(sl37);%向量ssl是針對信號s1構(gòu)造的向量disp=('正常信號的輸出向量')ssl＝[sl0,s11,sl2,sl3,s14,s15,sl6,s17]%＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝%用db1小波包對故障信號s2進行三層分解[t,d]=wpdec(s2,3,'db1','shannon');%plottree(t)%畫小波包樹結(jié)構(gòu)的圖形%s230是指信號S2的[3,0]結(jié)點的重構(gòu)系數(shù),其他以此類推s230=wprcoef(t,d,[3,0]);s231=wprcoef(t,d,[3,1]);s232=wprcoef(t,d,[3,2]);s233=wprcoef(t,d,[3,3]);s234=wprcoef(t,d,[3,4]);s235=wprcoef(t,d,[3,5]);s236=wprcoef(t,d,[3,6]);s237=wprcoef(t,d,[3,7]);%畫出重構(gòu)系數(shù)的波形subplot(9,2,4);plot(s230);Ylabel('S230');subplot(9,2,6);plot(s231);Ylabel('S231');subplot(9,2,8);plot(s232);Ylabel('S232');subplot(9,2,10);plot(s233);Ylabel('S233');subplot(9,2,12);plot(s234);Ylabel('S234');subplot(9,2,14);plot(s235);Ylabel('S235');subplot(9,2,16);plot(s236);Ylabel('S236');subplot(9,2,18);plot(s237);Ylabel('S237');%----------------------------------------------------------%計算故障信號各重構(gòu)系數(shù)的方差%s20是指s230的方差，其他依次類推s20＝norm(s230);s21＝norm(s231);s22＝norm(s232);s23＝norm(s233);s24＝norm(s234);s25＝norm(s235);s26＝norm(s236);s27＝norm(s237);%向量ss2是針對信號S1構(gòu)造的向量disp('故障信號')來源于：\o"軸承故障診斷小波分析.7z_小波軸承故障-專業(yè)指導(dǎo)其他資源-CSDN下載"軸承故障診斷小波分析.7z_小波軸承故障-專業(yè)指導(dǎo)其他資源-CSDN下載這問的方法實在太多了，代碼也有很全的，重點就是寫論文上。關(guān)于問題二：1、PCA故障檢測clc;clear;%%1.導(dǎo)入數(shù)據(jù)%產(chǎn)生訓(xùn)練數(shù)據(jù)num_sample=100;a=10*randn(num_sample,1);x1=a+randn(num_sample,1);x2=1*sin(a)+randn(num_sample,1);x3=5*cos(5*a)+randn(num_sample,1);x4=0.8*x2+0.1*x3+randn(num_sample,1);xx_train=[x1,x2,x3,x4];%產(chǎn)生測試數(shù)據(jù)a=10*randn(num_sample,1);x1=a+randn(num_sample,1);x2=1*sin(a)+randn(num_sample,1);x3=5*cos(5*a)+randn(num_sample,1);x4=0.8*x2+0.1*x3+randn(num_sample,1);xx_test=[x1,x2,x3,x4];xx_test(51:100,2)=xx_test(51:100,2)+15*ones(50,1);%%2.數(shù)據(jù)處理Xtrain=xx_train;Xtest=xx_test;X_mean=mean(Xtrain);X_std=std(Xtrain);[X_row,X_col]=size(Xtrain);Xtrain=(Xtrain-repmat(X_mean,X_row,1))./repmat(X_std,X_row,1);%標(biāo)準(zhǔn)化處理%%3.PCA降維SXtrain=cov(Xtrain);%求協(xié)方差矩陣[T,lm]=eig(SXtrain);%求特征值及特征向量,特征值排列順序為從小到大D=flipud(diag(lm));%將特征值從大到小排列%確定降維后的數(shù)量num=1;whilesum(D(1:num))/sum(D)<0.85num=num+1;endP=T(:,X_col-num+1:X_col);%取對應(yīng)的向量P_=fliplr(P);%特征向量由大到小排列%%4.計算T2和Q的限值%求置信度為99%時的T2統(tǒng)計控制限,T=k*(n^2-1)/n(n-k)*F(k,n-k)%其中k對應(yīng)num,n對應(yīng)X_rowT2UCL1=num*(X_row-1)*(X_row+1)*finv(0.99,num,X_row-num)/(X_row*(X_row-num));%求置信度為99%時的T2統(tǒng)計控制限%求置信度為99%的Q統(tǒng)計控制限fori=1:3th(i)=sum((D(num+1:X_col)).^i);endh0=1-2*th(1)*th(3)/(3*th(2)^2);ca=norminv(0.99,0,1);QU=th(1)*(h0*ca*sqrt(2*th(2))/th(1)+1+th(2)*h0*(h0-1)/th(1)^2)^(1/h0);%置信度為99%的Q統(tǒng)計控制限%%5.模型測試n=size(Xtest,1);Xtest=(Xtest-repmat(X_mean,n,1))./repmat(X_std,n,1);%標(biāo)準(zhǔn)化處理%求T2統(tǒng)計量，Q統(tǒng)計量[r,y]=size(P*P');I=eye(r,y);T2=zeros(n,1);Q=zeros(n,1);lm_=fliplr(flipud(lm));%T2的計算公式Xtest.T*P_*inv(S)*P_*Xtestfori=1:nT2(i)=Xtest(i,:)*P_*inv(lm_(1:num,1:num))*P_'*Xtest(i,:)';Q(i)=Xtest(i,:)*(I-P*P')*Xtest(i,:)';end%%6.繪制T2和SPE圖figure('Name','PCA');subplot(2,1,1);plot(1:i,T2(1:i),'k');holdon;plot(i:n,T2(i:n),'k');title('統(tǒng)計量變化圖');xlabel('采樣數(shù)');ylabel('T2');holdon;line([0,n],[T2UCL1,T2UCL1],'LineStyle','--','Color','r');subplot(2,1,2);plot(1:i,Q(1:i),'k');holdon;plot(i:n,Q(i:n),'k');title('統(tǒng)計量變化圖');xlabel('采樣數(shù)');ylabel('SPE');holdon;line([0,n],[QU,QU],'LineStyle','--','Color','r');%%7.繪制貢獻圖%7.1.確定造成失控狀態(tài)的得分S=Xtest(51,:)*P(:,1:num);r=[];fori=1:numifS(i)^2/lm_(i)>T2UCL1/numr=cat(2,r,i);endend%7.2.計算每個變量相對于上述失控得分的貢獻cont=zeros(length(r),X_col);fori=length(r)forj=1:X_colcont(i,j)=abs(S(i)/D(i)*P(j,i)*Xtest(51,j));endend%7.3.計算每個變量的總貢獻CONTJ=zeros(X_col,1);forj=1:X_colCONTJ(j)=sum(cont(:,j));end%7.4.計算每個變量對Q的貢獻e=Xtest(51,:)*(I-P*P');%選取第60個樣本來檢測哪個變量出現(xiàn)問題。contq=e.^2;%5.繪制貢獻圖figuresubplot(2,1,1);bar(contq,'g');xlabel('變量號');ylabel('SPE貢獻率%');holdon;subplot(2,1,2);bar(CONTJ,'r');xlabel('變量號');ylabel('T^2貢獻率%');來源于：\o"基于PCA的故障診斷方法（matlab）_湯憲宇的博客-CSDN博客_pca故障診斷"基于PCA的故障診斷方法（matlab）_湯憲宇的博客-CSDN博客_pca故障診斷不懂的轉(zhuǎn)原文看，花個時間學(xué)學(xué)。2、《普通公路隧道機電設(shè)施故障檢測方法研究》中的方法這個方法就是計算，看懂這篇論文就行。3、小波分析這個方法也看論文學(xué)習(xí)，找找代碼。代碼需要下載\o"使用小波進行信號中的特征檢測：使用小波進行特征檢測的MATLAB代碼和相關(guān)文件-matlab開發(fā)_-互聯(lián)網(wǎng)文檔類資源-CSDN下載"使用小波進行信號中的特征檢測：使用小波進行特征檢測的MATLAB代碼和相關(guān)文件-matlab開發(fā)_-互聯(lián)網(wǎng)文檔類資源-CSDN下載\o"基于小波分析的故障檢測與診斷碩士論文-基于小波分析的故障檢測與診斷.rar_-互聯(lián)網(wǎng)文檔類資源-CSDN下載"基于小波分析的故障檢測與診斷碩士論文-基于小波分析的故障檢測與診斷.rar_-互聯(lián)網(wǎng)文檔類資源-CSDN下載兩篇都有，我也沒下過，但看標(biāo)題來說。第一個是代碼，第二個是論文。至于性能評價無非就是那幾個指標(biāo)，準(zhǔn)確率、召回率等等。%樣本標(biāo)記為0和1，num為選取前n個特征的數(shù)據(jù)用于分類%需要安裝好SVMfunction[sens,spec,F1,pre,rec,acc]=SEERES(train,trainclass,test,testclass,num)acc=zeros(num,1);sens=zeros(num,1);spec=zeros(num,1);F1=zeros(num,1);pre=zeros(num,1);rec=zeros(num,1);FeatureNumber=zeros(num,1);[len,b]=size(testclass);forn=1:numlabel=trainclass;data=train(:,1:n);testlabel=testclass;testdata=test(:,1:n);model=svmtrain(label,data,'-s0-t0-b1');%默認C-SVC類型，0線性2RBF,-b會輸出概率[predictlabel,accuracy,Scores]=svmpredict(testlabel,testdata,model,'-b1');acc(n,1)=accuracy(1,1);FeatureNumber(n,1)=n;tp=0;fn=0;fp=0;tn=0;fory=1:lenifpredictlabel(y,1)==1&&testclass(y,1)==1tp=tp+1;elseifpredictlabel(y,1)==1&&testclass(y,1)==0fp=fp+1;elseifpredictlabel(y,1)==0&&testclass(y,1)==1fn=fn+1;elseifpredictlabel(y,1)==0&&testclass(y,1)==0tn=tn+1;endendsens(n,1)=tp/(tp+fn);spec(n,1)=tn/(tn+fp);pre(n,1)=tp/(tp+fp);rec(n,1)=sens(n,1);F1(n,1)=2*(pre(n,1)*rec(n,1))/(pre(n,1)+rec(n,1));end來源于：\o"查準(zhǔn)率、召回率、敏感性、特異性和F1-score的計算及Matlab實現(xiàn)_黑山白雪m的博客-CSDN博客_召回率和敏感性"查準(zhǔn)率、召回率、敏感性、特異性和F1-score的計算及Matlab實現(xiàn)_黑山白雪m的博客-CSDN博客_召回率和敏感性關(guān)于問題三：1、LSTM故障診斷\o"基于LSTM的故障診斷_舊日之歌的博客-CSDN博客_故障檢測lstm"基于LSTM的故障診斷_舊日之歌的博客-CSDN博客_故障檢測lstm看看上面的這個博文。2、小波分析與神經(jīng)網(wǎng)絡(luò)這個方法的代碼在CSDN里需要下載，但是這個很成熟了，而且已經(jīng)有人用這個方法做過診斷齒輪箱的故障。大家可以去看看下面這篇論文?！痘谛〔ń翟牒虰P神經(jīng)網(wǎng)絡(luò)的風(fēng)力發(fā)電機組齒輪箱故障診斷研究》就是代碼比較難找，但方法應(yīng)該是比較好的。3、深度學(xué)習(xí)故障診斷算法這個方法python厲害，擅長編程的可以考慮，用的是殘差收縮網(wǎng)絡(luò)。#!/usr/bin/envpython3#-*-coding:utf-8-*-"""CreatedonSatDec2823:24:052019ImplementedusingTensorFlow1.0.1andKeras2.2.1M.Zhao,S.Zhong,X.Fu,etal.,DeepResidualShrinkageNetworksforFaultDiagnosis,IEEETransactionsonIndustrialInformatics,2019,DOI:10.1109/TII.2019.2943898@author:super_9527"""from__future__importprint_functionimportkerasimportnumpyasnpfromkeras.datasetsimportmnistfromkeras.layersimportDense,Conv2D,BatchNormalization,Activationfromkeras.layersimportAveragePooling2D,Input,GlobalAveragePooling2Dfromkeras.optimizersimportAdamfromkeras.regularizersimportl2fromkerasimportbackendasKfromkeras.modelsimportModelfromkeras.layers.coreimportLambdaK.set_learning_phase(1)#Inputimagedimensionsimg_rows,img_cols=28,28#Thedata,splitbetweentrainandtestsets(x_train,y_train),(x_test,y_test)=mnist.load_data()ifK.image_data_format()=='channels_first':x_train=x_train.reshape(x_train.shape[0],1,img_rows,img_cols)x_test=x_test.reshape(x_test.shape[0],1,img_rows,img_cols)input_shape=(1,img_rows,img_cols)else:x_train=x_train.reshape(x_train.shape[0],img_rows,img_cols,1)x_test=x_test.reshape(x_test.shape[0],img_rows,img_cols,1)input_shape=(img_rows,img_cols,1)#Noiseddatax_train=x_train.astype('float32')/255.+0.5*np.random.random([x_train.shape[0],img_rows,img_cols,1])x_test=x_test.astype('float32')/255.+0.5*np.random.random([x_test.shape[0],img_rows,img_cols,1])print('x_trainshape:',x_train.shape)print(x_train.shape[0],'trainsamples')print(x_test.shape[0],'testsamples')#convertclassvectorstobinaryclassmatricesy_train=keras.utils.to_categorical(y_train,10)y_test=keras.utils.to_categorical(y_test,10)defabs_backend(inputs):returnK.abs(inputs)defexpand_dim_backend(inputs):returnK.expand_dims(K.expand_dims(inputs,1),1)defsign_backend(inputs):returnK.sign(inputs)defpad_backend(inputs,in_channels,out_channels):pad_dim=(out_channels-in_channels)//2inputs=K.expand_dims(inputs,-1)inputs=K.spatial_3d_padding(inputs,((0,0),(0,0),(pad_dim,pad_dim)),'channels_last')returnK.squeeze(inputs,-1)#ResidualShrinakgeBlockdefresidual_shrinkage_block(incoming,nb_blocks,out_channels,downsample=False,downsample_strides=2):residual=incomingin_channels=incoming.get_shape().as_list()[-1]foriinrange(nb_blocks):identity=residualifnotdownsample:downsample_strides=1residual=BatchNormalization()(residual)residual=Activation('relu')(residual)residual=Conv2D(out_channels,3,strides=(downsample_strides,downsample_strides),padding='same',kernel_initializer='he_normal',kernel_regularizer=l2(1e-4))(residual)residual=BatchNormalization()(residual)residual=Activation('relu')(residual)residual=Conv2D(out_channels,3,padding='same',kernel_initializer='he_normal',kernel_regularizer=l2(1e-4))(residual)#Calculateglobalmeansresidual_abs=Lambda(abs_backend)(residual)abs_mean=GlobalAveragePooling2D()(residual_abs)#Calculatescalingcoefficientsscales=Dense(out_channels,activation=None,kernel_initializer='he_normal',kernel_regularizer=l2(1e-4))(abs_mean)scales=BatchNormalization()(scales)scales=Activation('relu')(scales)scales=Dense(out_channels,activation='sigmoid',kernel_regularizer=l2(1e-4))(scales)scales=Lambda(expand_dim_backend)(scales)#Calculatethresholdsthres=keras.layers.multiply([abs_mean,scales])#Softthresholdingsub=keras.layers.subtract([residual_abs,thres])zeros=keras.layers.subtract([sub,sub])n_sub=keras.layers.maximum([sub,zeros])residual=keras.layers.multiply([Lambda(sign_backend)(residual),n_sub])#DownsamplingusingthepooL-sizeof(1,1)ifdownsample_strides>1:identity=AveragePooling2D(pool_size=(1,1),strides=(2,2))(identity)#Zero_paddingtomatchchannelsifin_channels!=out_channels:identity=Lambda(pad_backend,arguments={'in_channels':in_channels,'out_channels':out_channels})(identity)residual=keras.layers.add([residual,identity])returnresidual#defineandtrainamodelinputs=Input(shape=input_shape)net=Conv2D(8,3,padding='same',kernel_initializer='he_normal',kernel_regularizer=l2(1e-4))(inputs)net=residual_shrinkage_block(net,1,8,downsample=True)net=BatchNormalization()(net)net=Activation('relu')(net)net=GlobalAveragePooling2D()(net)outputs=Dense(10,activation='softmax',kernel_initializer='he_normal',kernel_regularizer=l2(1e-4))(net)model=Model(inputs=inputs,outputs=outputs)pile(loss='categorical_crossentropy',optimizer=Adam(),metrics=['accuracy'])model.fit(x_train,y_train,batch_size=100,epochs=5,verbose=1,validation_data=(x_test,y_test))#getresultsK.set_learning_phase(0)DRSN_train_score=model.evaluate(x_train,y_train,batch_size=100,verbose=0)print('Trainloss:',DRSN_train_score[0])print('Trainaccuracy:',DRSN_train_score[1])DRSN_test_score=model.evaluate(x_test,y_test,batch_size=100,verbose=0)print('Testloss:',DRSN_test_score[0])print('Testaccuracy:',DRSN_test_score[1])來源于：\o"深度學(xué)習(xí)故障診斷算法：殘差收縮網(wǎng)絡(luò)_weixin_47174159的博客-CSDN博客_故障診斷算法"深度學(xué)習(xí)故障診斷算法：殘差收縮網(wǎng)絡(luò)_weixin_47174159的博客-CSDN博客_故障診斷算法性能評價借鑒前一問的。關(guān)于第四問：直接運用前面的模型就行，注意論文寫好點。第三問代碼的更新：LSTM故障診斷matlab代碼：其中第一個用于讀取并劃分原始數(shù)據(jù)

第二個用于完成劃分訓(xùn)練集測試集，特征提取+分類等工作//本函數(shù)Input為//interval-數(shù)據(jù)劃分長度，默認為6400，即每6400個數(shù)據(jù)點劃為一個樣本//ind_column-通道，默認為2，即選取第二個通道//Output為//label1label2,...,label8，分別為劃分好的8種故障類型的樣本%function[label1label2label3label4label5label6label7label8]=read_data_1800_High(interval,ind_column)ifnargin<2ind_column=2;//如果傳遞的實參小于2個，默認ind_column為2endifnargin<1interval=6400;//默認interval=6400endfile_rul='E:\Datasets\PHMdatachallenge\2009PHMSocietyConferenceDataChallenge-gearbox\spur_30hz_High\';//以下為獲取file_rul路徑下.mat格式的所有文件file_folder=fullfile(file_rul);dir_output=dir(fullfile(file_folder,'*.mat'));file_name={dir_}';num_file=max(size(file_name));//num_file為文件數(shù)，本例中num_file=8，8個文件，分別存儲齒輪箱的8種故障數(shù)據(jù)fori=1:num_filefile=[file_rul,file_name{i}];load(file);[filepath,name,ext]=fileparts(file);raw=eval(name); //每6400個點劃分為一個樣本n=1;left_index=1+(n-1)*interval;right_index=n*interval;whileright_index<=size(raw,1)temp=raw(left_index:right_index,ind_column);//eval函數(shù)構(gòu)造label1,label2,...等變量名eval(['label'num2str(i)'(:,n)=temp;']);n=n+1;left_index=1+(n-1)*interval;right_index=n*interval;endendend//讀取數(shù)據(jù)，label1為一個6400*83的數(shù)組，83為每種故障類型所得到的樣本數(shù)[label1,label2,label3,label4,label5,label6,label7,label8]=read_data_1800_High();num_categories=8;//由于matlab中LSTM建模需要，用num2cell函數(shù)將label1轉(zhuǎn)為cell型，label_x_cell為一個1×83的cell型數(shù)組，每個cell存儲6400個數(shù)據(jù)點label1_x_cell=num2cell(label1,1);label2_x_cell=num2cell(label2,1);label3_x_cell=num2cell(label3,1);label4_x_cell=num2cell(label4,1);label5_x_cell=num2cell(label5,1);label6_x_cell=num2cell(label6,1);label7_x_cell=num2cell(label7,1);label8_x_cell=num2cell(label8,1);num_1=length(label1_x_cell);num_2=length(label2_x_cell);num_3=length(label3_x_cell);num_4=length(label4_x_cell);num_5=length(label5_x_cell);num_6=length(label6_x_cell);num_7=length(label7_x_cell);num_8=length(label8_x_cell);//創(chuàng)建用于存儲每種故障類型的標(biāo)簽的數(shù)據(jù)結(jié)構(gòu)，由于matlab中l(wèi)stm建模需要，也需要cell型數(shù)據(jù)。例如，label1_y為一個83×1的cell型數(shù)組，目前其值為空label1_y=cell(num_1,1);label2_y=cell(num_2,1);label3_y=cell(num_3,1);label4_y=cell(num_4,1);label5_y=cell(num_5,1);label6_y=cell(num_6,1);label7_y=cell(num_7,1);label8_y=cell(num_8,1);//創(chuàng)建故障類型的標(biāo)簽，用1,2,3,...,8表示8種故障標(biāo)簽，給對應(yīng)標(biāo)簽賦值。fori=1:num_1;label1_y{i}='1';endfori=1:num_2;label2_y{i}='2';endfori=1:num_3;label3_y{i}='3';endfori=1:num_4;label4_y{i}='4';endfori=1:num_5;label5_y{i}='5';endfori=1:num_6;label6_y{i}='6';endfori=1:num_7;label7_y{i}='7';endfori=1:num_8;label8_y{i}='8';end//用dividerand函數(shù)將每種故障類型的數(shù)據(jù)隨機劃分為4:1的比例，分別用作訓(xùn)練和測試[trainInd_label1,~,testInd_label1]=dividerand(num_1,0.8,0,0.2);[trainInd_label2,~,testInd_label2]=dividerand(num_2,0.8,0,0.2);[trainInd_label3,~,testInd_label3]=dividerand(num_3,0.8,0,0.2);[trainInd_label4,~,testInd_label4]=dividerand(num_4,0.8,0,0.2);[trainInd_label5,~,testInd_label5]=dividerand(num_5,0.8,0,0.2);[trainInd_label6,~,testInd_label6]=dividerand(num_6,0.8,0,0.2);[trainInd_label7,~,testInd_label7]=dividerand(num_7,0.8,0,0.2);[trainInd_label8,~,testInd_label8]=dividerand(num_8,0.8,0,0.2);//構(gòu)建每種故障類型的訓(xùn)練數(shù)據(jù)xTrain_label1=label1_x_cell(trainInd_label1);yTrain_label1=label1_y(trainInd_label1);xTrain_label2=label2_x_cell(trainInd_label2);yTrain_label2=label2_y(trainInd_label2);xTrain_label3=label3_x_cell(trainInd_label3);yTrain_label3=label3_y(trainInd_label3);xTrain_label4=label4_x_cell(trainInd_label4);yTrain_label4=label4_y(trainInd_label4);xTrain_label5=label5_x_cell(trainInd_label5);yTrain_label5=label5_y(trainInd_label5);xTrain_label6=label6_x_cell(trainInd_label6);yTrain_label6=label6_y(trainInd_label6);xTrain_label7=label7_x_cell(trainInd_label7);yTrain_label7=label7_y(trainInd_label7);xTrain_label8=label8_x_cell(trainInd_label8);yTrain_label8=label8_y(trainInd_label8);//構(gòu)建每種故障類型的測試數(shù)據(jù)xTest_label1=label1_x_cell(testInd_label1);yTest_label1=label1_y(testInd_label1);xTest_label2=label2_x_cell(testInd_label2);yTest_label2=label2_y(testInd_label2);xTest_label3=label3_x_cell(testInd_label3);yTest_label3=label3_y(testInd_label3);xTest_label4=label4_x_cell(testInd_label4);yTest_label4=label4_y(testInd_label4);xTest_label5=label5_x_cell(testInd_label5);yTest_label5=label5_y(testInd_label5);xTest_label6=label6_x_cell(testInd_label6);yTest_label6=label6_y(testInd_label6);xTest_label7=label7_x_cell(testInd_label7);yTest_label7=label7_y(testInd_label7);xTest_label8=label8_x_cell(testInd_label8);yTest_label8=label8_y(testInd_label8);//將每種故障類型的數(shù)據(jù)整合，構(gòu)建完整的訓(xùn)練集和測試集xTrain=[xTrain_label1xTrain_label2xTrain_label3xTrain_label4xTrain_label5xTrain_label6xTrain_label7xTrain_label8];yTrain=[yTrain_label1;yTrain_label2;yTrain_label3;yTrain_label4;yTrain_label5;yTrain_label6;yTrain_label7;yTrain_label8];num_train=size(xTrain,2);xTest=[xTest_label1xTest_label2xTest_label3xTest_label4xTest_label5xTest_label6xTest_label7xTest_label8];yTest=[yTest_label1;yTest_label2;yTest_label3;yTest_label4;yTest_label5;yTest_label6;yTest_label7;yTest_label8];num_test=size(xTest,2);//================================================================================//以下分別對每個樣本，提取三種特征：1.瞬時頻率，2.瞬時譜熵，3.小波包能量，//上述三種特征后面會被送入分類器中進行分類，實驗結(jié)果表明，將小波包能量作為特征，//能夠取得最高的分類精度//提取瞬時頻率：用matlab的pspectrum對每個樣本進行譜分解，再用instfreq函數(shù)計算瞬時頻率FreqResolu=25;TimeResolu=0.12;//theoutputofpspectrum'p'containsanestimateoftheshort-term,time-localizedpowerspectrumofx.//Inthiscase,pisofsizeNf×Nt,whereNfisthelengthoffandNtisthelengthoft.[p,f,t]=cellfun(@(x)pspectrum(x,fs,'TimeResolution',TimeResolu,'spectrogram'),xTrain,'UniformOutput',false);instfreqTrain=cellfun(@(x,y,z)instfreq(x,y,z)',p,f,t,'UniformOutput',false);[p,f,t]=cellfun(@(x)pspectrum(x,fs,'TimeResolution',TimeResolu,'spectrogram'),xTest,'UniformOutput',false);instfreqTest=cellfun(@(x,y,z)instfreq(x,y,z)',p,f,t,'UniformOutput',false);//提取瞬時譜熵：用matlab的pspectrum對每個樣本進行譜分解，再用pentropy函數(shù)計算瞬時頻率[p,f,t]=cellfun(@(x)pspectrum(x,fs,'TimeResolution',TimeResolu,'spectrogram'),xTrain,'UniformOutput',false);pentropyTrain=cellfun(@(x,y,z)pentropy(x,y,z)',p,f,t,'UniformOutput',false);[p,f,t]=cellfun(@(x)pspectrum(x,fs,'TimeResolution',TimeResolu,'spectrogram'),xTest,'UniformOutput',false);pentropyTest=cellfun(@(x,y,z)pentropy(x,y,z)',p,f,t,'UniformOutput',false);//提取小波包能量//num_level=5表示進行小波包五層分解，共獲得2^5=32個值組成的特征向量。num_level=5;index=0:1:2^num_level-1;//wpdec為小波包分解函數(shù)treeTrain=cellfun(@(x)wpdec(x,num_level,'dmey'),xTrain,'UniformOutput',false);treeTest=cellfun(@(x)wpdec(x,num_level,'dmey'),xTest,'UniformOutput',false);fori=1:num_trainforj=1:length(index) //wprcoef為小波系數(shù)重構(gòu)函數(shù)reconstr_coef=wprcoef(treeTrain{i},[num_level,index(j)]);//計算能量energy(j)=sum(reconstr_coef.^2);endenergyTrain_doule(i,:)=energy;endenergyTrain=num2cell(energyTrain_doule,2);energyTrain=energyTrain';fori=1:num_testforj=1:length