matlab_时域统计特征_特征提取_轴承_

%% Script to test the Empirical Wavelet Transform % Based on the papers: % J. Gilles, "Empirical Wavelet Transform", IEEE % Trans. on Signal Processing, 2013 % J. Gilles, G. Tran, S. Osher, "2D Empirical transforms. Wavelets, % Ridgelets and Curvelets Revisited", SIAM Journal on Imaging Sciences, % 2014 % J. Gilles, K. Heal, "A parameterless scale-space approach to find % meaningful modes in histograms - Application to image and spectrum % segmentation", submitted 2014. % % Don't hesitate to modify the parameters and try % your own signals! % % Author: Jerome Gilles % Institution: UCLA - Department of Mathematics % Year: 2014 % Version: 2.0 clear all %% User setup load('E:\毕业设计\毕设\轴承数据\Normal Baseline Data\98.mat') a3=X098_DE_time(1:61440,1); b3=reshape(a3,2048,30)'; signal = b3(1,:)'; f=signal; t=0:length(f)-1; % Choose the signal you want to analyze % (sig1,sig2,sig3,sig4=ECG,sig5=seismic,sig6=EEG) % signal = 'sig4'; params.SamplingRate = 12000; %put -1 if you don't know the sampling rate %params.SamplingRate = 4000; %put -1 if you don't know the sampling rate % channel = 50; %for EEG only % Choose the wanted global trend removal (none,plaw,poly,morpho,tophat) params.globtrend = 'none'; params.degree=6; % degree for the polynomial interpolation 与选择poly有关 % Choose the wanted regularization (none,gaussian,avaerage,closing) params.reg = 'none'; params.lengthFilter = 10; %与选择gaussian有关，高斯掩模的长度 params.sigmaFilter = 1.5; %与选择gaussian有关，高斯标准偏差 % Choose the wanted detection method (locmax,locmaxmin,ftc, % adaptive,adaptivereg,scalespace) params.detect ='locmaxmin'; % params.detect = 'scalespace'; % params.typeDetect='otsu'; %for scalespace:otsu,halfnormal,empiricallaw,mean,kmeans params.N = 5; % maximum number of bands params.completion = 0; % choose if you want to force to have params.N modes % in case the algorithm found less ones (0 or 1) %params.InitBounds = [4 8 13 30]; % params.InitBounds = [2 25]; %与选择adaptive,adaptivereg有关 % Perform the detection on the log spectrum instead the spectrum params.log=0; % Choose the results you want to display (Show=1, Not Show=0) Bound=1; % Display the detected boundaries on the spectrum Comp=1; % Display the EWT components Rec=0; % Display the reconstructed signal TFplane=0; % Display the time-frequency plane (by using the Hilbert % transform). You can decrease the frequency resolution by % changing the subresf variable below. Demd=1; % Display the Hilbert-Huang transform (YOU NEED TO HAVE % FLANDRIN'S EMD TOOLBOX) subresf=1; % InitBounds = params.InitBounds; % switch lower(signal) % case 'sig1' % load('sig1.mat'); % t=0:1/length(f):1-1/length(f); % case 'sig2' % load('sig2.mat'); % t=0:1/length(f):1-1/length(f); % case 'sig3' % load('sig3.mat'); % t=0:1/length(f):1-1/length(f); % case 'sig4' % load('sig4.mat'); % t=0:length(f)-1; % case 'sig5' % load('seismic.mat'); % f=f(10000:20000); %sub portion of the signal used in the paper % t=0:length(f)-1; % case 'sig6' % load('eeg.mat') % t=0:length(f)-1; % end %% We perform the empirical transform and its inverse % compute the EWT (and get the corresponding filter bank and list of % boundaries) [ewt,mfb,boundaries]=EWT1D(f,params); %% Show the results if Bound==1 %Show the boundaries on the spectrum div=1; if (strcmp(params.detect,'adaptive')||strcmp(params.detect,'adaptivereg')) Show_EWT_Boundaries(abs(fft(f)),boundaries,div,params.SamplingRate,InitBounds); else Show_EWT_Boundaries(abs(fft(f)),boundaries,div,params.SamplingRate); end end if Comp==1 %Show the EWT components and the reconstructed signal if Rec==1 %compute the reconstruction rec=iEWT1D(ewt,mfb); Show_EWT(ewt,f,rec); else Show_EWT(ewt); end end % if TFplane==1 %Show the time-frequency plane by using the Hilbert transform % Hilbert_EWT(ewt,t,f,1,0); % end if TFplane==1 %Show the time-frequency plane by using the Hilbert transform EWT_TF_Plan(ewt,boundaries,params.SamplingRate,f,[],[],subresf,[]); end for i=1:5 %b=kurt(ewt{i,:},'kurt2'); b(i)=kurtosis(ewt{i,:}); % %b=kurtosis(imf(i,:)); % a=ewt{i,:}; % %a=imf(i,:); g = corrcoef(f,a); g1(i)=g(1,2); % n=2; % ystd=std(a,1); % r=0.2*ystd; % h= approx_entropy_my(n,r,a); % h end Data2=ewt{5,:}+ewt{3,:}; fs=12000; N=2048; % ts = 0:1/fs:(N-1)/fs; % 采样时刻 fre=fs*(0:N/2-1)/N; %真实频率，(0:N/2-1)/N为前一半的序列,可以看到，频率的最大值为fs/2，频率的分辨率为fs/Ns=1/T。 y=fft(Data2,N); mag_y=abs(y)/(N/2); %这里N取得是信号的有效长度，如果在信号末尾补0的话，不算做有效长度。 figure plot(fre,mag_y(1:N/2)); % xlim([0 1000]) xlabel('频率/Hz'); ylabel('振幅'); title('信号的频谱'); %% EMD comparison: if you have Patrick Flandrin's EMD toolbox you can %% perform the EMD and display the corresponding Time-Frequency plane if Demd==1 imf=emd(f); Disp_HHT(imf,t,f,1,1); end