figure(3) subplot(5,1,1) plot(imf(1,:)); xlabel('时间/min'); ylabel('IMF1/MW');
时间: 2023-11-18 14:04:50 浏览: 23
这段代码的作用是在一个包含5个子图的图形中,绘制经过 EMD 处理得到的第一分量(IMF1)随时间的变化趋势。其中,subplot(5,1,1) 表示将当前图形分割成5行1列,当前绘制的是第1个子图。plot(imf(1,:)) 表示绘制 imf 中第1行的数据,也就是第一分量随时间的变化趋势。xlabel('时间/min') 和 ylabel('IMF1/MW') 则是为绘图添加横轴和纵轴的标签。
相关问题
优化一下下面这段MATLAB代码:clear; clc; S=load('ExSignal.txt'); fs=2000; t=0:(1/fs):0.4; [m,n]=size(S); signal=S(1:m,n); % 对信号进行频谱分析 N = 2^nextpow2(fs); Y = fft(signal, N)/N*2; f = fs/N*(0:1:N-1); A = abs(Y); P = angle(Y); figure; subplot(211); plot(t, signal); title('滤波前信号'); xlabel('时间(s)'); ylabel('幅值'); subplot(212); plot(f(1:N/2), A(1:N/2)); title('幅值频谱'); xlabel('频率(Hz)'); ylabel('幅值'); % 设计带通滤波器1 d1 = fdesign.bandpass('Fst1,Fp1,Fp2,Fst2,Ast1,Ap,Ast2', 30, 35, 60, 65, 60, 1, 60, fs); h1 = design(d1); % 对信号进行带通滤波并画出波形 newsignal = filter(h1, signal); Y1 = fft(newsignal, N)/N*2; A1 = abs(Y1); figure; subplot(211); plot(t, newsignal); title('带通滤波后信号'); xlabel('时间(s)'); ylabel('幅值'); subplot(212); plot(f(1:N/2), A1(1:N/2)); title('幅值频谱'); xlabel('频率(Hz)'); ylabel('幅值'); % 设计带通滤波器2 d2 = fdesign.bandpass('Fst1,Fp1,Fp2,Fst2,Ast1,Ap,Ast2', 230, 235, 285, 290, 60, 1, 60, fs); h2 = design(d2); % 对信号进行带通滤波并画出波形 newsignal = filter(h2, signal); Y1 = fft(newsignal, N)/N*2; A1 = abs(Y1); figure; subplot(211); plot(t, newsignal); title('带通滤波后信号'); xlabel('时间(s)'); ylabel('幅值'); subplot(212); plot(f(1:N/2), A1(1:N/2)); title('幅值频谱'); xlabel('频率(Hz)'); ylabel('幅值');
优化后的代码如下:
clear;
clc;
% 加载信号
S = load('ExSignal.txt');
signal = S(:, end); % 取最后一列
% 设定采样率和时间轴
fs = 2000;
t = 0:1/fs:0.4;
% 对信号进行频谱分析
N = 2^nextpow2(length(signal));
Y = fft(signal, N)/N*2;
f = fs/N*(0:1:N-1);
A = abs(Y);
% 绘制原始信号和幅值频谱
figure;
subplot(211);
plot(t, signal);
title('滤波前信号');
xlabel('时间(s)');
ylabel('幅值');
subplot(212);
plot(f(1:N/2), A(1:N/2));
title('幅值频谱');
xlabel('频率(Hz)');
ylabel('幅值');
% 设计带通滤波器1
d1 = fdesign.bandpass('Fst1,Fp1,Fp2,Fst2,Ast1,Ap,Ast2', 30, 35, 60, 65, 60, 1, 60, fs);
h1 = design(d1);
% 对信号进行带通滤波并画出波形和幅值频谱
newsignal1 = filter(h1, signal);
Y1 = fft(newsignal1, N)/N*2;
A1 = abs(Y1);
figure;
subplot(211);
plot(t, newsignal1);
title('带通滤波1后信号');
xlabel('时间(s)');
ylabel('幅值');
subplot(212);
plot(f(1:N/2), A1(1:N/2));
title('幅值频谱');
xlabel('频率(Hz)');
ylabel('幅值');
% 设计带通滤波器2
d2 = fdesign.bandpass('Fst1,Fp1,Fp2,Fst2,Ast1,Ap,Ast2', 230, 235, 285, 290, 60, 1, 60, fs);
h2 = design(d2);
% 对信号进行带通滤波并画出波形和幅值频谱
newsignal2 = filter(h2, signal);
Y2 = fft(newsignal2, N)/N*2;
A2 = abs(Y2);
figure;
subplot(211);
plot(t, newsignal2);
title('带通滤波2后信号');
xlabel('时间(s)');
ylabel('幅值');
subplot(212);
plot(f(1:N/2), A2(1:N/2));
title('幅值频谱');
xlabel('频率(Hz)');
ylabel('幅值');
% 清理无用变量
clearvars -except signal newsignal1 newsignal2;
clear %清除内存 load('1797 b007_0.mat') %根据实际需要更改地址、路径 sig=X118_DE_time(1:12000); fs=12000; N=12000; Ts=1/fs; sig=sig(1:N);%设置取样频率fs,取样数N t=0:Ts:(N-1)*Ts;%时间轴 t sig=(sig-mean(sig))/std(sig,1);%对 sig 进行归一化 subplot(211);plot(t,sig);%绘制 sig 波形 xlabel('时间 t/s'); ylabel('振动加速度/V'); nfft=fs/2; % 16384 S=pspectrum(sig,nfft);%对 sig 做功率谱 subplot(212); plot((0:nfft/2 -1)/nfft*fs,S(1:nfft/2));% 绘制功率谱 xlabel('频率 f/Hz'); ylabel('功率谱 P/W') [c,l]=wavedec(sig,3,'db2');%利用 db2 对 sig 进行 3 级小波分解 c3=wrcoef ('a',c ,l,'db2',3); d3=wrcoef('d',c,l,'db2',3); d2 =wrcoef ('d',c,l,'db2',2); d1 =wrcoef('d',c,l,'db2',1);%重构第 1-3 层细节 d1~d3 和第 3 层概貌 c3 figure; subplot(414); plot(t,c3); ylabel('c3');%绘制 c3 subplot(413); plot(t,d3); ylabel('d3');%绘制 d3 subplot(412); plot(t,d2); ylabel('d2');%绘制 d2 subplot(411); plot(t,d1); ylabel('d1');%绘制 d1 y=hilbert(d1); %对 d1 进行 Hilbert 变换,得y ydata=abs(y); %ydata=|y| ydata=ydata-mean(ydata);%对 ydata 去均值(目的是去除幅度较大的直流分量) P=pspectrum(ydata,nfft);%ydata 的功率谱为 P figure; plot((0:nfft/2-1)/nfft*fs,P(1:nfft/2)); xlabel('频率 f/Hz');%绘出 d1 的 Hilbert 包络谱 P=P(1:nfft/2); [M,f1]=max(P); f1=f1*fs/nfft-1 %故障频率 f1为包络谱中幅度最大处的频率 将代码由利用db2进行3级小波分解改为利用db10进行5级小波分解
clear %清除内存
load('1797 b007_0.mat') %根据实际需要更改地址、路径
sig=X118_DE_time(1:12000);
fs=12000;
N=12000;
Ts=1/fs;
sig=sig(1:N);%设置取样频率fs,取样数N
t=0:Ts:(N-1)*Ts;%时间轴
t
sig=(sig-mean(sig))/std(sig,1);%对 sig 进行归一化
subplot(211);
plot(t,sig);%绘制 sig 波形
xlabel('时间 t/s');
ylabel('振动加速度/V');
nfft=fs/2; % 16384
S=pspectrum(sig,nfft);%对 sig 做功率谱
subplot(212);
plot((0:nfft/2 -1)/nfft*fs,S(1:nfft/2));% 绘制功率谱
xlabel('频率 f/Hz');
ylabel('功率谱 P/W')
[c,l]=wavedec(sig,5,'db10');%利用 db10 对 sig 进行 5 级小波分解
c5=wrcoef ('a',c ,l,'db10',5);
d5=wrcoef('d',c,l,'db10',5);
d4 =wrcoef ('d',c,l,'db10',4);
d3 =wrcoef('d',c,l,'db10',3);
d2 =wrcoef ('d',c,l,'db10',2);
d1 =wrcoef('d',c,l,'db10',1);%重构第 1-5 层细节 d1~d5 和第 5 层概貌 c5
figure;
subplot(511); plot(t,sig); ylabel('原始信号');%绘制 sig
subplot(512); plot(t,c5); ylabel('c5');%绘制 c5
subplot(513); plot(t,d5); ylabel('d5');%绘制 d5
subplot(514); plot(t,d4); ylabel('d4');%绘制 d4
subplot(515); plot(t,d3); ylabel('d3');%绘制 d3
y=hilbert(d1); %对 d1 进行 Hilbert 变换,得y
ydata=abs(y); %ydata=|y|
ydata=ydata-mean(ydata);%对 ydata 去均值(目的是去除幅度较大的直流分量)
P=pspectrum(ydata,nfft);%ydata 的功率谱为 P
figure;
plot((0:nfft/2-1)/nfft*fs,P(1:nfft/2));
xlabel('频率 f/Hz');%绘出 d1 的 Hilbert 包络谱
P=P(1:nfft/2);
[M,f1]=max(P);
f1=f1*fs/nfft-1 %故障频率 f1为包络谱中幅度最大处的频率