信号处理-MATLAB-共振解调-CWRU数据集

%共振解调,轴承故障分析程序。Bearing failure analysis %以下随便选取CWRU四个文本，并结合说明文件中的Defect frequencies计算理论值(f=Df*rpm/60),具体如下 % 97.mat正常数据文本 0HP(Horse Power，就是负载) % 105.mat 内圈故障 7mil 0HP 159.92Hz % 198.mat 外圈故障 14mil 1HP 105.87Hz % 3005.mat 滚动体故障 28mil 0HP 141.7Hz % 另外文件夹中还有3006-3008另外三个滚动体故障文本，X049_DE_time-X051_DE_time，可更换验证共振解调算法 % 3006 1HP 1797rpm 139.2Hz % 3007 2HP 1772rpm 137.4Hz % 3008 3HP 1730rpm 135.9Hz clc; clear; close all; %参数设置 sampleFreq = 12000;%采样频率 sampleLength = 8196;%采样点数 %下载正常数据，数据爱放哪放哪，反正MATLAB能识别中文路径，路径别过长就行 good = load('D:\\\97.mat'); data_good = good.X097_DE_time(1:sampleLength); %下载内圈故障数据 inner = load('D:\\\105.mat'); data_inner = inner.X105_DE_time(1:sampleLength); %下载外圈故障数据 outer = load('D:\\\198.mat'); data_outer = outer.X198_DE_time(1:sampleLength); %下载滚动体故障数据 ball = load('D:\\\3005.mat'); data_ball = ball.X048_DE_time(1:sampleLength); %产生仿真信号 t=0:1/sampleFreq:(sampleLength-1)/sampleFreq; %信号赋值 Data1 = data_ball; %========================================================================== %直接进行FFT变换 %========================================================================== for n=1:sampleLength if abs(Data1(n))>10; Data(n)=0; else Data(n)=Data1(n); end end fft_result = abs(fft(Data)) * 2 / sampleLength; %画图的坐标变换 time_plot_s = 0:1/sampleFreq:(sampleLength-1) / sampleFreq; fft_plot_Hz = sampleFreq*(1:sampleLength/2)/sampleLength; figure; subplot(211) plot(time_plot_s,Data); title('时域波形');ylabel('振幅/m/s^2');xlabel('时间/s'); grid on; subplot(212) plot(fft_plot_Hz,fft_result(1:sampleLength/2)); title('频域波形');ylabel('振幅/m/s^2');xlabel('频率/Hz'); grid on; %当butter(4,[0.01,0.1])，将下边预处理两部分注释掉或者保留都可以，加不加预处理都不会对共振解调结果产生明显影响 %========================================================================== %预处理：去除趋势项 %========================================================================== m=1; data=Data'; a=polyfit(t,Data,m); y=Data-polyval(a,t); %========================================================================== %预处理：小波降噪 %========================================================================== [c,l]=wavedec(y,3,'db3'); [thr,sorh,keepapp]=ddencmp('den','wv',y); Data=wdencmp('gbl',y,'db3',3,thr,sorh,keepapp); %========================================================================== %共振解调法 %========================================================================== %带通滤波 subplot(211); plot(time_plot_s,Data); % 包络 envelop_hil = hilbert(Data); envelop_abs = abs(envelop_hil); [b,a] = butter(4,[0.01,0.1]) ; filter_data = filter(b,a,envelop_abs); %fft变换 envelop_fft = abs(fft(filter_data))*2 /sampleLength; figure; subplot(211) plot(time_plot_s,envelop_abs); title('包络时域波形');ylabel('振幅/m/s^2');xlabel('时间/s');grid on; subplot(212) plot(fft_plot_Hz,10*envelop_fft(2:sampleLength/2+1)); title('包络频域波形');ylabel('振幅/m/s^2');xlabel('频率/Hz');grid on; % N=length(envelop_fft); % Xrms=sqrt(sum(envelop_fft.^2)/N) % beta=sum(envelop_fft.^4)/N % kv=beta/(Xrms^4)