EMD算法的matlab程序介绍解析

%此版本为ALAN版本的整合注释版

fun ctio n imf =emd(x

%Empiricial Mode Decompositi on (Hilbert-Hua ngTra nsform

%imf =emd(x

%Func :fin dpeaks

x=tra nspose(x(:;%转置为行矩阵

imf =[];

while ~ismonotonic(x%当x不是单调函数，分解终止条件

x1=x;

sd =lnf;% 均值

%直到x1满足IMF条件得c1

while (sd>0.1 |~isimf(x1%当标准偏差系数sd大于0.1或x1不是固有模态函数时，分量终止条件

s仁getspli ne(x1;%上包络线

s2=-getspli ne(-x1;% 下包络线

x2=x1-(s1+s2/2;%此处的x2为文章中的h

sd =sum((x1-x2.A2/sum(x1.A2;

x1=x2;

end

imf{e nd+1}=x1;

x =x-x1;

end

imf{en d+1}=x;

%FUNCTIONS

fun ctio n u =ism onotonic(x

%u=0表示x不是单调函数,u=1表示x为单调的

u1=le ngth(fi ndpeaks(x*le ngth(fi ndpeaks(-x;

if u1>0, u =0;

else, u =1; end

function u =isimf(x

%u=0表示x不是固有模式函数，u=1表示x是固有模式函数N =le ngth(x;

u仁sum(x(1:N-1.*x(2:N<0;

u2=le ngth(fi ndpeaks(x+le ngth(fi ndpeaks(-x;

if abs(u1-u2>1, u =0;

else, u =1; end

function s =getspli ne(x

%三次样条函数拟合成元数据包络线

N =le ngth(x;

p =fin dpeaks(x;

s =spline([0p N+1],[Ox(pO],1:N;

function n =fin dpeaks(x

%Fi nd peaks.找到极值,n为极值点所在位置%n =fin dpeaks(x

n =fin d(diff(diff(x>0 <0;

u =fin d(x (n+1>x( n;

n(u=n( u+1;

function plot_hht00(x,Ts

%双边带调幅信号的EMD分解

%Plot the HHT.

%plot_hht(x,Ts

%

%::Sy ntax

%The array 列x is the in put sig nal and Ts is the sampli ng period 取样周期.%Example on use:[x,Fs]=wavread('Hum.wav';

%plot_hht(x(1:6000,1/Fs;

%Func :emd

%Get HHT.

clear all;

close all;

Ts=0.0005;

t=0:Ts:10;% 采样率2000HZ

%调幅信号

%x=si n(2*pi*t.*si n(40*pi*t;

x=si n(2*pi*t;

s仁getspl in e(x;%上包络线

s2=-getspli ne(-x;% 上包络线

x仁(s1+s2/2;%此处的x2为文章中的h

figure;

plot(t,x;xlabel('Time',ylabel('Amplitude';title('双边带调幅信号';hold on;

plot(t,s1,'-r';

plot(t,s2,'-r';

y 轴句柄

plot(t,x1,'g';

imf =emd(x;

for k =1:le ngth(imf

b(k=sum(imf{k}.*imf{k};

th =an gle(hilbert(imf{k};

d{k}=diff(th/Ts/(2*pi;

end

[u,v]=sort(-b;

b =1-b/max(b;

%Set time-freque ncy plots.

N =le ngth(x;

c =li nspace(0,(N-2*Ts,N-1;

%

figure;

for k =v(1:2

plot(c,d{k},'k.','Color',b([kk k],'MarkerSize',3;hold on; set(gca,'FontSize',8,'XLim',[0c(end],'YLim',[050];% 设置 x 、 xlabel('Time',ylabel('Frequency';title('原信号时频图’；

end

%Set IMF plots.

M =le ngth(imf;

N =le ngth(x;

c =li nspace(0,(N-1*Ts,N;

for k1=0:4:M-1

figure

for k2=1:m in (4,M-k1,

subplot(4,1,k2,

plot(c,imf{k1+k2};

set(gca,'Fo ntSize',8,'XLim',[0c(e nd]; title('EMD分解结果';

end

xlabel('Time';

end