轨道交通装备健康管理五次作业

Homework_1 Data Preprocessing --Wavelet Transform题目：Original raw signal is provided. Choose an appropriated wavelet function with levels. Then conduct signal preprocessing using wavelet de -noising. Post your Matlab code and plot a result figure such as:Matlab 程序：clear all clcs = load('sta_ch1_1.txt'); % 加载原始信号dn = length(s); % length of vector or largest array dimension[C,L] = wavedec(s,6,'db5');% Help file, C: the wavelet decomposition vector , L: the bookkeeping vectorA6 = wrcoef('a', C, L,'db5',6); % a: approximation (基本信号)低频部分（尺度系数，逼近系数）D1 = wrcoef('d', C, L, 'db5',1); % d: detail (细节信号)高频部分（小波系数，细节系数）D2 = wrcoef('d', C, L, 'db5',2); D3 = wrcoef('d', C, L, 'db5',3); D4 = wrcoef('d', C, L, 'db5',4); D5 = wrcoef('d', C, L, 'db5',5);D6 = wrcoef('d', C, L, 'db5',6);subplot(811); plot(s,'k','LineWidth',1.5); ylabel('s','FontSize',12,'Fontweight','bold'); title('sta_ch1_1');subplot(812); plot(A6(1:dn/2^1),'k','LineWidth',1.5); ylabel('a6','FontSize',12,'Fontweight','bold')subplot(813); plot(D1(1:dn/2^1),'k','LineWidth',1.5);ssta c h11a 1d 1d2d 3d 4d5d 6ylabel('d1','FontSize',12,'Fontweight','bold')subplot(814); plot(D2(1:dn/2^2),'k','LineWidth',1.5); ylabel('d2','FontSize',12,'Fontweight','bold')subplot(815); plot(D3(1:dn/2^3),'k','LineWidth',1.5); ylabel('d3','FontSize',12,'Fontweight','bold')subplot(816); plot(D4(1:dn/2^4),'k','LineWidth',1.5); ylabel('d4','FontSize',12,'Fontweight','bold')subplot(817); plot(D5(1:dn/2^5),'k','LineWidth',1.5); ylabel('d5','FontSize',12,'Fontweight','bold')subplot(818); plot(D6(1:dn/2^6),'k','LineWidth',1.5); ylabel('d6','FontSize',12,'Fontweight','bold')运行结果：Homework_2 Batch data processing题目：1) Matlab code program to calculate at least 3 features from time -domain and frequency -domain respectively;2) Project calculated features on figures of 2-D and 3-D with good discrimination looking like:2-D Plot with one selected feature 3-D Plot with three selected features Key: Batch data processing and saveMatlab 程序：clear all clcna=dir('*_first_*.txt'); % lists the files N=size(na,1);% plot3(m(:,1),m(:,2),m(:,3),'b.') for k=1:Ns(:,k) = load(na(k).name); % 加载原始信号 end% s=decimalscaling(s);% [s,ps]=mapminmax(s',0.1,.9); % s=s';% s=mapminmax('reverse',s',ps); % s=s';F1F 3Number of dataN o .7 f e a t u r emymean=mean(s);mystd=std(s);myskewness=skewness(s);mykurtosis=kurtosis(s);myrms=rms(s);shapefactor=myrms./mean(s);crestfactor=max(s)./myrms;figure(1);subplot(321);plot(mymean,'b.');xlabel('mean');subplot(322);plot(mystd,'b.');xlabel('std');subplot(323);plot(myskewness,'b.');xlabel('skewness');subplot(324);plot(mykurtosis,'b.');xlabel('kurtosis');subplot(325);plot(shapefactor,'b.');xlabel('shapefactor');subplot(326);plot(crestfactor,'b.');xlabel('crestfactor');% figure(2);plot3(mymean,mystd,myskewness,'b.')% figure(3);plot3(mymean,mystd,mykurtosis,'b.')% figure(4);plot3(mymean,myskewness,mykurtosis,'b.')figure(5);plot3(mystd,mymean,shapefactor,'b.');xlabel('std');ylabel(' mean');zlabel('shapefactor');运行结果：Homework_3 Feature Extraction by PCA题目：A data matrix with 30 samples in rows and 168 features in colums1) Matlab programing to conduct features dimension reduction by PCA method, post your codes on assignment.2) Choose three of new transformed features and project them onto a 2D or 3D space shown like:Hints:1) To select the dimensions number k of new features, setupImportance of Principal ComponentsPrincipal ComponentsE n g e n v a l u eNew F1Features project after PCANew F2N e w F 3cumulative 0.8 of variance captured by k dimensions.No of PC’s1234567∑∑===mi iki i NoPCs 11λλ2) Suggested matlab functions: svd, repmat, find, sort, and princomp.Matlab 程序：clc clearX = load('train_Flux2.dat'); [row,col]=size(X); mymean=mean(X);X2=X-repmat(mymean,row,1);%减去平均值归一化[U,S,V] = svd(X2,0);%U basis vectors S singular values V eigen values totalPCs=sum(diag(S)); temp=0; k=0;rate=0; while rate<0.8 k=k+1;temp=temp+S(k,k);rate=temp/totalPCs; endlength(find(S~=0)); figure(1)bar(S(S~=0));hold on %S(find(S~=0))title('Contribution of Principal Components')xlabel('Principal Components');ylabel('PCs importance'); y=0:0.01:S(1);plot([k,k],[0,S(1)],'-r','linewidth',2)figure(2)V=V(1:end,1:k);X3=X*V;%新的特征30*10矩阵，30组数据每组都有十个新特征 plot3(X3(:,1),X3(:,2),X(:,5),'.'); 运行结果：Homework_4 Sequential Probability Ratio Test(SPRT)题目：Given a simulated time -domain signal with injected anomaly after 900 points (Fig. 1), make Matlab programing to realize SPRT anomaly detection with H2 hypothesis test, and plot the result like Fig.2.SetupSimulated time-domain signalTime StateS P R T I n d e x1)Cut data points 101 to 700 as healthy dataset, cut data points701 to 1000 as test dataset;2)false alarm probabilityα=0.001; miss alarm probabilityβ=0.005,Mean shift M = standard deviation of healthy dataset;Matlab程序：clearclose allclcData = load('data.dat');Test = Data(701:1000,:); % 101 to 700 healthy,701 to 1000 test dataset; t = size(Test,1); % 测试样本数alpha=0.001; %虚警率beta=0.005; %漏警率A=log(beta./(1-alpha)); %判断的下限值，低于为正常B=log((1-beta)./alpha); %判断的上限值，超过为故障mean_normal=mean(Data(101:700));sigma=std(Data(101:700),1);mu_abn=mean(Test(201:300));i=1;SPRT_index(1)=(mu_abn-mean_normal)*Test(1)/sigma^2-(mu_abn^2-mean_nor mal^2)/(2*sigma^2);while i<tif SPRT_index(i)<B && SPRT_index(i)>Ai= i+1;SPRT_index(i)=SPRT_index(i-1)+(mu_abn-mean_normal)*Test(i)/sigma^2-(mu_abn^2-mean_n ormal^2)/(2*sigma^2);elsei =i+1;SPRT_index(i) =(mu_abn-mean_normal)*Test(i)/sigma^2-(mu_abn^2-mean_normal^2)/(2*sigm a^2);endendfor i=1:tif SPRT_index(i)>=BSPRT_index(i)=B; %似然比对数值大于上限值，则赋值为上限值Belseif SPRT_index(i)<=ASPRT_index(i)=A; %似然比对数值小于下限值，则赋值为下限值Aend%介于之间，则不变if SPRT_index(i)>0SPRT_index(i)=SPRT_index(i)./B; %若值为正，则除以Belse%否则，除以A的绝对值SPRT_index(i)= SPRT_index(i)./abs(A);endendfigure(1)plot(SPRT_index,'*--black'); hold onalarm_label = find(SPRT_index >= 1);alarm_warning = SPRT_index(alarm_label);plot(alarm_label,alarm_warning,'ro','LineWidth',2)title('Fault detection - H2 hypothesis');xlabel('Time State')ylabel('SPRT Index')运行结果：。