emd程序.txt

谐波小波程序变参数的积分程序多点人工生成地震波程序低通滤波程序经验模式分解降噪程序

N = 1024;% # of data samples
fs=1000
t=(1:N)/fs;
%t=1:1024
%x=sin(2*pi*8*t)+sin(2*pi*20*t)
%x=cos(2*pi*4*t+0.5*sin(2*pi*2*t))+sin(2*pi*16*t);
x=(1+0.5*sin(2*pi*1*t)).*cos(2*pi*14*t+0.4*sin(2*pi*3*t))+sin(2*pi*2*t)+sin(2*pi*7*t);
imf=emd(x);
%[A,fa,tt]=hhspectrum(imf);
if size(imf,1)>1
    [A,fa,tt]=hhspectrum(imf(1:end-1,:));;
else
    [A,fa,tt]=hhspectrum(imf);
end
[E,tt1]=toimage(A,fa,tt,length(tt));
for k=1:size(E,1)
    bjp(k)=sum(E(k,:))*1/fs;
end
f=(0:N-3)/N*(fs/2);
%%绘制原信号
figure(1)
plot(t,x);
xlabel('时间');
ylabel('幅值')
%%%%%%EMD分解
figure(2)
[m,n]=size(imf);
for i=1:m    
    subplot(m,1,i);
    plot(imf(i,:));
    axis([0,1000,-2,2])
    xlabel('时间');
    ylabel('频率')
end
%emd_visu(x,t,imf,1);
%%%%绘制灰度图
figure(3)
imagesc(tt1,[0,0.5*fs],E);
set(gca,'YDir','normal')
colormap(flipud(gray))
axis([0,1000,0,50])
xlabel('频率 / Hz');
ylabel('幅值')
%%%绘制边际谱
figure(4)
plot(f,2*bjp)
axis([0,50,0,1.2])
xlabel('频率 / Hz');
ylabel('幅值')
%绘制mesh三维图,不可用
figure(5)
mesh(tt,fa,abs(A))
xlabel('time/s')
ylabel('frequency/Hz')
zlabel('Ampititude')
title('[WV] mesh')