clear;clc;close all;
Ts = (1/2048);
fs = 1/Ts;
t=0:Ts:1-Ts;
freqs = 2*pi*(t-0.5-1/Ts)/(fs);
os1=sin(2*pi*30*t+pi/2); %高频正弦信号
os2=sin(2*pi*8*t+pi/3); %低频正弦信号
os3=[zeros(1,300),1.5*rand(1,600),zeros(1,300),1.5*rand(1,500),zeros(1,348)];%间歇信号
os=os1+os2+os3;
x=os;
% x=awgn(os,10,'measured','dB');
T=length(x);
%绘制仿真信号和其频谱图
figure
subplot(411);plot(t,os1, 'r');xlabel('t/s');ylabel('幅值');legend('信号1');
title('sin(60*pi*t)')
hold on
subplot(412);plot(t,os2, 'r');xlabel('t/s');ylabel('幅值');legend('信号2');
title('sin(16*pi*t)')
hold on
subplot(413);
plot(t,os3, 'r');xlabel('t/s');ylabel('幅值');legend('间歇信号');
hold on
subplot(414);
plot(t,x, 'r');xlabel('t/s');ylabel('幅值');legend('含噪间歇信号');
figure
y2=x;
L=length(y2);
NFFT = 2^nextpow2(L);
Y = fft(y2,NFFT)/L;
f = fs/2*linspace(0,1,NFFT/2);
plot(f,2*abs(Y(1:NFFT/2)));
xlabel('f/Hz');ylabel('|PF(f)|');legend('含噪间歇信号的频谱');
% LMD分解
[pf,a,si,u] = lmd(x);
line=size(pf,1)
NN = length(pf(1,:))
t = linspace(0,1,NN);
for k1=0:4:line-1
figure('Color',[1 1 1]);
for k2=1:min(4,line-k1)
subplot(4,2,2*k2-1);
plot(t,pf(k1+k2,:));
title(sprintf('第%d个PF', k1+k2))
xlabel('Time/s')
ylabel(sprintf('PF%d',k1+k2));legend('LMD');
subplot(4,2,2*k2)
[yf, f] = FFTAnalysis(pf(k1+k2,:), Ts);
plot(f, yf)
title(sprintf('第%d个PF的频谱', k1+k2))
xlabel('f/Hz')
ylabel('|PF(f)|');legend('LMD');
end
end;
for i = 1:1:line
cc(i)=min(min(corrcoef(pf(i,:), x)))
end
figure
plot(cc,'-g<','LineWidth',1.5,'MarkerEdgeColor','b','MarkerFaceColor','b','MarkerSize',5);
set(gca,'XGrid', 'on', 'YGrid', 'on');
legend('CC');
xlabel('IMF');
ylabel('LMD相关系数');
%重构信号
cg_ev1=pf(1,:)+pf(2,:)+pf(3,:);
figure;
plot(t,x,'-b');hold on
plot(t,cg_ev1,'-r');xlabel('t/s');ylabel('幅值');legend('原始信号','LMD重构信号');
%误差信号
err1=os-cg_ev1;
figure
plot(t,err1);xlabel('t/s');ylabel('幅值');legend('LMD误差信号');
% ELMD分析
Nstd=0.2;
NR=200;
modes = mlmd(x,Nstd,NR);
PF=modes;
line=size(PF,1);
NN = length(PF(1,:));
n = linspace(0,1,NN);
for k1=0:5:line-1
figure('Color',[1 1 1]);
for k2=1:min(5,line-k1)
subplot(5,2,2*k2-1);
plot(t,PF(k1+k2,:));
title(sprintf('第%d个PF', k1+k2))
xlabel('Time/s')
ylabel(sprintf('PF%d',k1+k2));legend('ELMD');
subplot(5,2,2*k2)
[yf, f] = FFTAnalysis(PF(k1+k2,:), Ts);
plot(f, yf)
title(sprintf('第%d个PF的频谱', k1+k2))
xlabel('f/Hz')
ylabel('|PF(f)|');legend('ELMD');
end
end;
for i = 1:1:line
cc(i)=min(min(corrcoef(modes(i,:), x)));
end
figure
plot(cc,'-g<','LineWidth',1.5,'MarkerEdgeColor','b','MarkerFaceColor','b','MarkerSize',5);
set(gca,'XGrid', 'on', 'YGrid', 'on');
legend('CC');
xlabel('IMF');
ylabel('ELMD的相关系数');
%重构信号
cg_ev2=PF(3,:)+PF(4,:);
figure;
subplot(211);plot(t,os);xlabel('t/s');ylabel('幅值');legend('原始信号');
hold on
subplot(212);
plot(t,cg_ev2);xlabel('t/s');ylabel('幅值');legend('ELMD的重构信号');
%误差信号
err2=os-cg_ev2;
figure
plot(t,err2);xlabel('t/s');ylabel('幅值');legend('ELMD误差信号');
%性能评价
Ps_two=sum(sum((os).^2));%signal power
Pn_two=sum(sum((cg_ev1-os).^2));%noise power
SRN1=10*log10(Ps_two/Pn_two)
MSE1=(sum(sum((os-cg_ev1).^2))/length(x)).^0.5
MAE1= mean(abs(os-cg_ev1))
CC1= min(min(corrcoef(os, cg_ev1)))
Ps_two=sum(sum((os).^2));%signal power
Pn_two=sum(sum((os-cg_ev2).^2));%noise power
SRN2=10*log10(Ps_two/Pn_two)
MSE2=(sum(sum((os-cg_ev2).^2))/length(x)).^0.5
MAE2= mean(abs(os-cg_ev2))
CC2= min(min(corrcoef(os, cg_ev2)))
- 1
- 2
- 3
前往页