上一个简单的EMD程序

空杯子

因为本人仍为新手，这学期才开始研究HHT，找了很多文献和程序，觉得下面的这个比较简单，适合初级者使用。已经加了一些解释，希望能更易懂。


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


plot_hht00(x,Ts)只画出了时频图，没有三维图。

[ 本帖最后由 空杯子 于 2009-5-15 16:19 编辑 ]

空杯子

本来想传仿真图的，更形象一些，但是无奈不会。。。

空杯子

jiangyuan881018

我最近做毕业论文，刚刚才设计这个领域的东西，你传的EMD对我很有用。谢谢了～不过，我对MATLAB实在是不熟练。我有了一些文本数据和这个程序后，我应该怎么让MATLAB实现啊？从而得到你下面的这些图啊？谢谢你了～

jiangyuan881018

另外，程序里的X代表是什么啊？是文本数据么？我是菜鸟，你可能觉得我问的问题很幼稚，可是我真的不懂！

空杯子

x代表你要EMD分解的信号，根据要求自己设定的

jyy

我最近也在弄毕业设计呢，我的课题是基于HHT的非平稳信号分析
我们老师要我编程仿真，我完全不懂

[ 本帖最后由 zhangnan3509 于 2009-8-18 21:20 编辑 ]

wenchengcq

请教：我的matlab软件没有程序中的getspline()函数，该怎么办？

wenchengcq

运行plot_hht00(x,Ts)时，程序一直在画图，不停地显示，

yizimeng

请问一下楼主这个程序中需要哪些数据,以及整个程序的主程序是怎样的啊??我正在学matlb,希望能够接触一个完整的程序,谢谢楼主了!!

sunsun1212

我最近做毕业论文，刚刚才设计这个领域的东西，你传的EMD对我很有用。谢谢你了

花月落

不错啊，非常好用啊

花月落

不错啊，非常好用

花月落

不错啊，非常好用

花月落

不错啊，非常好用