|
楼主 |
发表于 2016-10-14 16:46
|
显示全部楼层
clear;clc;clf;
N=1000;
%fft默认计算的信号是从0开始的
Ts=1/500;
t=(1:N)*Ts;
x=sin(2*pi*3*t)+0.4*sin(2*pi*11*t)+0.9*sin(2*pi*25*t);
z=x;
c=emd(z);
%计算每个IMF分量及最后一个剩余分量residual与原始信号的相关性
[m,n]=size(c);
for i=1:m;
a=corrcoef(c(i,:),z);
xg(i)=a(1,2);
end
xg;
for i=1:m-1
%--------------------------------------------------------------------
%计算各IMF的方差贡献率
%定义:方差为平方的均值减去均值的平方
%均值的平方
%imfp2=mean(c(i,:),2).^2
%平方的均值
%imf2p=mean(c(i,:).^2,2)
%各个IMF的方差
mse(i)=mean(c(i,:).^2,2)-mean(c(i,:),2).^2;
end;
mmse=sum(mse);
for i=1:m-1
mse(i)=mean(c(i,:).^2,2)-mean(c(i,:),2).^2;
%方差百分比,也就是方差贡献率
mseb(i)=mse(i)/mmse*100;
%显示各个IMF的方差和贡献率
end;
%画出每个IMF分量及最后一个剩余分量residual的图形
figure(1)
for i=1:m-1
disp(['imf',int2str(i)]) ;disp([mse(i) mseb(i)]);
end;
subplot(m+1,1,1)
plot(z)
set(gca,'fontname','times New Roman')
set(gca,'fontsize',14.0)
ylabel(['signal','Amplitude'])
for i=1:m-1
subplot(m+1,1,i+1);
set(gcf,'color','w')
plot(c(i,:),'k')
set(gca,'fontname','times New Roman')
set(gca,'fontsize',14.0)
ylabel(['imf',int2str(i)])
end
subplot(m+1,1,m+1);
set(gcf,'color','w')
plot(c(m,:),'k')
set(gca,'fontname','times New Roman')
set(gca,'fontsize',14.0)
ylabel(['r',int2str(m-1)])
for i = 1:m-1
MiTest(i)=mi(z',c(i,:)');
end
MaxMI=max(MiTest);
for i = 1:m-1
NormalMI(i)=MiTest(i)/MaxMI;
end
|
|