emd 希尔伯特黄变换程序

（一）简单的EMD程序

function imf = emd(x)

% Empiricial Mode Decomposition (Hilbert-Huang Transform)

% imf = emd(x)

% Func : findpeaks

x = transpose(x(:));%转置

imf = [];

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

x1 = x;

sd = Inf;%均值

%直到x1满足IMF条件，得c1

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

s1 = getspline(x1);%上包络线

s2 = -getspline(-x1);%下包络线

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

sd = sum((x1-x2).^2)/sum(x1.^2);

x1 = x2;

end

imf{end+1} = x1;

x = x-x1;

end

imf{end+1} = x;

% FUNCTIONS

function u = ismonotonic(x)

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

u1 = length(findpeaks(x))*length(findpeaks(-x));

if u1 > 0, u = 0;

else, u = 1; end

function u = isimf(x)

%u=0表示x不是固有模式函数，u=1表示x是固有模式函数

N = length(x);

u1 = sum(x(1:N-1).*x(2:N) < 0);

u2 = length(findpeaks(x))+length(findpeaks(-x));

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

else, u = 1; end

function s = getspline(x)

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

N = length(x);

p = findpeaks(x);

s = spline([0 p N+1],[0 x(p) 0],1:N);

-------------------------------------------------------------------------------- function n = findpeaks(x)

% Find peaks.找到极值

% n = findpeaks(x)

n = find(diff(diff(x) > 0) < 0);

u = find(x(n+1) > x(n));

n(u) = n(u)+1;

------------------------------------------------------------------------------------------ ---------------------------------------------------------------------------------------- function plot_hht00(x,Ts)

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

% Plot the HHT.

% plot_hht(x,Ts)

%

% :: Syntax

% The array x is the input signal and Ts is the sampling 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:1; % 采样率2000HZ

% 调幅信号

x=sin(2*pi*t).*sin(40*pi*t);

s1 = getspline(x);%上包络线

s2 = -getspline(-x);%上包络线

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

figure;

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

plot(t,s1,'-r');

plot(t,s2,'-r');

plot(t,x1,'g');

imf = emd(x);

for k = 1:length(imf)

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

th = angle(hilbert(imf{k}));

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

end

[u,v] = sort(-b);

b = 1-b/max(b);

% Set time-frequency plots.

N = length(x);

c = linspace(0,(N-2)*Ts,N-1);

%

figure;

for k = v(1:2)

plot(c,d{k},'k.','Color',b([k k k]),'MarkerSize',3); hold on;

set(gca,'FontSize',8,'XLim',[0 c(end)],'YLim',[0 50]);%设置x、y轴句柄

xlabel('Time'), ylabel('Frequency');title('原信号时频图');

end

% Set IMF plots.

M = length(imf);

N = length(x);

c = linspace(0,(N-1)*Ts,N);

for k1 = 0:4:M-1

figure

for k2 = 1:min(4,M-k1), subplot(4,1,k2), plot(c,imf{k1+k2});

set(gca,'FontSize',8,'XLim',[0 c(end)]);title('EMD分解结果');end

xlabel('Time');

end

(二)

function imf = emd(x)

% Empiricial Mode Decomposition (Hilbert-Huang Transform)

% imf = emd(x)

% Func : findpeaks

x = transpose(x(:));%转置

imf = [];

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

x1 = x;

sd = Inf;%均值

%直到x1满足IMF条件，得c1

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

s1 = getspline(x1);%上包络线

s2 = -getspline(-x1);%下包络线

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

sd = sum((x1-x2).^2)/sum(x1.^2);

x1 = x2;

end

imf{end+1} = x1;