程序1:原始信号时域分析及小波去噪处理
clear all
z=importdata('C:\\Users\\wangkun\\Desktop\\轴承诊断\\test2.mat'); x1=z.X105_DE_time(1:4096); clear z; N=4096; fs=12000; n=0:N-1; t=n/fs; f=n*fs/N; figure(1); plot(t,x1);
xlabel('t'); ylabel('幅值'); title('原信号时域图') %小波去噪
[thr,sorh,keepapp]=ddencmp('den','wv',x1); xd=wdencmp('gbl',x1,'db3',2,thr,sorh,keepapp); figure(2); plot(t,xd); xlabel('t');
ylabel('幅值');
title('小波去噪后时域图')
程序2:EMD分解及Hilbert包络
clc
clear all
z=load('C:\\Users\\wangkun\\Desktop\\轴承诊断\\test2.mat'); x=z.X105_DE_time(1:1024); N=1024; fs=12000; n=0:N-1; f=n*fs/N; lag=N; n=0:N-1;
t=n/fs; imf=emd(x);
[m,n]=size(imf); %imf为一m*n阶矩阵,m是imf分量,n为数据点
emd_visu(x,1:length(x),imf,m); %实信号的信号重构及emd结果显示函数 for i=1:m
a(i)=kurtosis(imf(i,:));%峭度 b(i)=mean(imf(i,:)); %均值;
c(i)=var(imf(i,:)); %方差; d(i)=std(imf(i,:)); %均方值
e(i)=std(imf(i,:)).^0.5; %均方根值 f(i)=skewness(imf(i,:)); %计算偏度 end
[k,c]=max(a); %k为峭度最大值,c为最大元素在数组中的位置 [r,lags]=xcorr(x,lag,'unbiased'); %计算序列的自相关函数 for i=1:m
[R,lags]=xcorr(imf(i,:),lag,'unbiased'); %计算序列的自相关函数
a=corrcoef(R(1:N/2),r(1:N/2)); %相关系数矩阵【对称】,主对角元素为1 xg(i)=abs(a(1,2)); %相关系数 end
[R,C]=max(xg); %R为最大值,C为最大元素在数组中的位置
figure(4);
y = hilbert(imf(C,:)); a = abs(y);%包络 b=fft(a); mag1=abs(b); mag=mag1*2/N; f1=(0:N-1)*fs/N;
plot(f1(1:N/2),mag(1:N/2)); %set(gca,'xlim',[0,.400]); title('包络'); xlabel('频率'); ylabel('幅值');
相关推荐: