EMD

G-Rilling EMD工具箱

http://blog.sina.com.cn/s/blog_6163bdeb0102e2cd.html《转载》

工具箱的安装

运行install_emd.m文件可以实现此工具箱的安装，uninstall_emd.m实现卸载。

 

安装中的问题

    但是安装的时候，如果使用的是VS的编译器（mbuild –setup、mex –setup设置），会报找不到complex.h的问题（用Linux下的Matlab不会出错），从而使cemdc2_fix.c等文件编译失败，这几个文件是为了快速实现计算EMD而用c编写的，所以即使编译失败，也不影响直接使用emd.m实现EMD功能。如果想编译成功，可如下修改：（摘自http://www.chinavib.com/thread-79866-1-1.html）

    G. Rilling 07年3月份的程序，运行作者的install_emd.m，出现找不到complex.h的问题，以下是个人的理解和解决过程：（个人的运行环境为matlab6.5）

complex.h的问题

    产生原因：采用matlab的C编译函数mex时，定义了C99_OK的宏（EMDS/make_emdc.m (28行)），利用的是ANSI C99标准如果个人的电脑中没有相关的支持，就会出现这个问题。

解决方法：EMDS/make_emdc.m中第28行中mex(’-DC99_OK‘,args(:))语句中的 '-DC99_OK' 即可。

    注意：

    改完之后，运行install_emd，会出现M_PI没有定义的问题，缺少了常数PI的宏定义，导致一些.c文件编译失败。

    产生原因：去掉C99_OK之后，程序中使用的是作者提供的 emd_complex.h和emd_complex.c两个文件来支持复数运算，这两个文件中，并没有定义M_PI这个宏。

    解决方法：M_PI这个宏，只在两个文件中（clocal_mean.c和clocal_mean2.c）使用，个人的解决方法是，在相应的头文件（clocal_mean.h和clocal_mean2.h）中加入M_PI的宏定义即可。

    在两个.h文件中分别加入一下语句：

#define CLOCAL_MEAN_H

#ifndef M_PI

#define M_PI 3.1415926

#endif

    安装完成后，编译输出的.dll文件会出现，重复后缀名的问题，及 xxx.dll 变成了 xxx.dll.dll自己去掉多余的.dll即可

    最后，关于版本问题：作者推荐使用7.1+版本，但只是针对个别的函数有影响，主要是作者提供的例子程序，无法在matlab6.5环境中运行，算法的主要功能函数并不受影响。

 

    如果懒得向上面修改，可以直接下载已经编译成功的工具箱压缩包，即里面已经有了cemdc.mexw32、cemdc_fix.mexw32等文件，这里可以下载到

http://www.chinavib.com/thread-84036-1-1.html

    下载后，通过设置路径，作为matlab工具箱使用即可。

 

工具箱的使用

工具箱函数

    运行help index_emd可以查看工具箱提供的函数，如下

index_emd.M list of functions in the EMD package
 
  type help function_name for more information on a specific function
 
  Empirical Mode Decomposition
 
    emd          - computes EMD and bivariate/complex EMD with various options
    emd_local    - computes local EMD variation
    emd_online   - computes on-line EMD variation. Note that it does not truly
                   apply on-line: the function is only a demonstration.
    emdc         - fast implementation for EMD with Cauchy-like stopping criterion
                   (requires compilation, see make_emdc function)
    emdc_fix     - fast implementation for EMD with predefined number of iterations
                   (requires compilation, see make_emdc function)
    cemdc        - fast implementation for bivariate/complex EMD (first algorithm)
                   with Cauchy-like stopping criterion (requires compilation,
                   see make_emdc function)
    cemdc_fix    - fast implementation for bivariate/complex EMD (first algorithm)
                   with predefined number of iterations (requires compilation,
                   see make_emdc function)
    cemdc2       - fast implementation for bivariate/complex EMD (second algorithm)
                   with Cauchy-like stopping criterion (requires compilation,
                   see make_emdc function)
    cemdc2_fix   - fast implementation for bivariate/complex EMD (second algorithm)
                   with predefined number of iterations (requires compilation,
                   see make_emdc function)
 
  Utilities
 
    install_emd   - setup Matlab's path and compile the C codes.
    uninstall_emd - revert the modifications made by install_emd and remove the
                    files (optional).
    make_emdc     - compile all C codes
    emd_visu      - visualization of EMD
    cemd_visu     - visualization of bivariate/complex EMD (automatically called
                    by emd_visu when the input is complex)
    cenvelope     - compute envelope curves for bivariate/complex EMD
    cemd_disp     - visualization of envelope curves and tube envelope
    plot3c        - plot a complex vector in 3 dimensions
    plotc         - plot the projection of a complex vector on a variable direction
    dirstretch    - directional stretching of a complex vector
    hhspectrum    - compute Hilbert-Huang spectrum (need the Time-Frequency Toolbox
                    http://tftb.nongnu.org)
    toimage       - transform a spectrum made of 1D functions (e.g., output of
                    "hhspectrum") in an 2D image
    disp_hhs      - display the image output of "toimage" as a Hilbert-Huang spectrum
    addtag        - add a tag to a graphic object (uses the Tag property as a list
                    of keywords or "tags")
    rmtag         - remove a tag from a graphic object (uses the Tag property as
                    a list of keywords or "tags")
    hastag        - test whether a graphic object has a specific tag (uses the Tag
                    property as a list of keywords or "tags")
    findtag       - find objects having a specific tag (uses the Tag property as
                    a list of keywords or "tags")
   
  Examples from G. Rilling, P. Flandrin and P. Gon鏰lves,
    "On Empirical Mode Decomposition and its algorithms"
    IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing
    NSIP-03, Grado (I), June 2003
   
    emd_fmsin         - Fig. 1: a 3-component example (need the Time-Frequency
                        Toolbox http://tftb.nongnu.org)
    emd_triang        - Fig. 2: another 3-component example
    emd_sampling      - Fig. 3: effect of sampling on 1 tone
    emd_separation    - Fig. 4: separation of 2 tones
    ex_online         - Sect 3.4: the way emd_online.m works
    triangular_signal - subroutine called by emd_triang (formerly triang.m)
   
  Examples from G. Rilling, P. Flandrin, P. Gon鏰lves and J. M. Lilly,
    "Bivariate Empirical Mode Decomposition",
    Signal Processing Letters (submitted)
 
    bivariate_EMD_principle        - Fig. 1: principle of the bivariate/complex EMD
    bivariate_EMD_mean_definitions - Fig. 2: definition of the mean for each algorithm.
                                     Also allows to test other signals and parameter sets.
    bivariate_EMD_illustration     - Fig. 3: illustration of the bivariate EMD
                                     on an oceanographic float position record

    稍做整理如下：

EMD分解函数

函数	功能
emd	计算EMD、双变量/复数EMD
emd_local	计算local EMD
emd_online	计算在线EMD（不是真正在线应用，此函数只是一个示范）
emdc	使用Cauchy-like停止准则的快速EMD实现，需编译
emdc_fix	使用预定义迭代次数的快速EMD实现，需编译
cemdc	使用Cauchy-like停止准则的快速双变量/复数EMD实现（方法1），需编译
cemdc_fix	使用预定义迭代次数的快速双变量/复数EMD实现（方法1），需编译
cemdc2	使用Cauchy-like停止准则的快速双变量/复数EMD实现（方法2），需编译
cemdc2_fix	使用预定义迭代次数的快速双变量/复数EMD实现（方法2），需编译

工具函数

函数	功能
install_emd	设置Matlab路径，编译c代码
uninstall_emd	回复install_emd做的修改，移除文件
make_emdc	编译c代码
emd_visu	EMD可视化
cemd_visu	双变量/复数EMD可视化（emd_visu的输入是双变量或复数时自动改为调用cemd_visu）
cenvelope	计算双变量EMD的包络曲线
cemd_disp	显示复数包络曲线
plot3c	三维中绘制复数向量
plotc	绘制复数向量在一个可变方向上的投影
dirstretch	复数向量的方向拉伸
hhspectrum	计算Hilbert-Huang谱（需要时频工具箱http://tftb.nongnu.org）
toimage	将一个一维函数谱转化为图像
disp_hhs	以Hilbert-Huang谱的形式显示toimage函数的输出
addtag	添加标签到一个图形对象
rmtag	移除标签从一个图形对象
hastag	测试一个图形对象是否有指定的标签
findtag	找有指定标签的图形对象

来自《On Empirical Mode Decomposition and its algorithms》的Examples

函数	功能
emd_fmsin	一个包含3组分的例子（需要时频工具箱）
emd_triang	另一个包含3组分的例子
emd_sampling	effect of sampling on 1 tone
emd_separation	separation of 2 tones
ex_online	the way emd_online.m works
triangular_signal	emd_triang文件调用的子程序

来自《Bivariate Empirical Mode Decomposition》的Examples

函数	功能
bivariate_EMD_principle	双变量/复数EMD原则
bivariate_EMD_mean_definitions	各种方法的平均值的定义
bivariate_EMD_illustration	双变量EMD在海洋漂流位置的应用图解

工具箱使用示例

EMD

clc

clear all

close all

% 原始数据

fs = 1000;

ts = 1/fs;

t=0:ts:0.3;

z=2*sin(2*pi*10*t) + 5.*sin(2*pi*100*t);

figure

plot(t, z)

title('原始信号')

% EMD

imf=emd(z);

emd_visu(z,t,imf)

[A,f,tt]=hhspectrum(imf);

[im,tt]=toimage(A,f);

disp_hhs(im);

边际谱

clc

clear all

close all

% 原始数据

fs = 1000;

ts = 1/fs;

t=0:ts:0.3;

y=2*sin(2*pi*10*t);% + 5.*sin(2*pi*100*t);

figure

plot(t, y)

title('原始信号')

% 求Hilbert-Huang谱

[A,fh,th] = hhspectrum(y);

figure

subplot(211)

plot(th*ts, A)

title('瞬时幅值') % 就是包络

subplot(212)

plot(th*ts, fh*fs)

title('瞬时频率')

% 显示结果

[im,tt,ff] = toimage(A,fh,th);

disp_hhs(im,tt)

colormap(flipud(gray))

% 编程实现显示

figure

imagesc(tt*ts,[0,0.5*fs],im);

ylabel('frequency/Hz')

set(gca,'YDir','normal')

xlabel('time/s')

title('Hilbert-Huang spectrum')

例子程序

   更详细的使用说明可以参见例子程序，如emd_fmsin.m程序，运行结果如下

EMD分解如下

   可以看到，EMD实现的3个组分的分离（即分别分解到了IMF1~3中），可见EMD的强大功能。

   为便于以后查阅，代码也贴过来吧：

N = 2000;% # of data samples
T = 1:4:N;
t = 1:N;

p = N/2;% period of the 2 sinusoidal FM's

% sinusoidal FM 1
fmin1 = 1/64;% min frequency
fmax1 = 1.5*1/8;% max frequency
x1 = fmsin(N,fmin1,fmax1,p,N/2,fmax1);

% sinusoidal FM 1
fmin2 = 1/32;% min frequency
fmax2 = 1.5*1/4;% max frequency
x2 = fmsin(N,fmin2,fmax2,p,N/2,fmax2);

% logon
f0 = 1.5*1/16;% center frequency
x3 = amgauss(N,N/2,N/8).*fmconst(N,f0);

a1 = 1;
a2 = 1;
a3 = 1;

x = real(a1*x1+a2*x2+a3*x3);
x = x/max(abs(x));

[imf,ort,nbits] = emd(x);

emd_visu(x,t,imf,1);

figure(1)

% time-frequency distributions
Nf = 256;% # of frequency bins
Nh = 127;% short-time window length
w = tftb_window(Nh,'Kaiser');

[s,rs] = tfrrsp(x,T,Nf,w,1);
[s,rs1] = tfrrsp(imf(1,:)',T,Nf,w,1);
[s,rs2] = tfrrsp(imf(2,:)',T,Nf,w,1);
[s,rs3] = tfrrsp(imf(3,:)',T,Nf,w,1);

figure(4)

subplot(221)
imagesc(flipud(rs(1:128,:)))
set(gca,'YTick',[]);set(gca,'XTick',[])
xlabel('time')
ylabel('frequency')
title('signal')

pause

subplot(222)
imagesc(flipud(rs1(1:128,:)))
set(gca,'YTick',[]);set(gca,'XTick',[])
xlabel('time')
ylabel('frequency')
title('mode #1')

pause

subplot(223)
imagesc(flipud(rs2(1:128,:)))
set(gca,'YTick',[]);set(gca,'XTick',[])
xlabel('time')
ylabel('frequency')
title('mode #2')

pause

subplot(224)
imagesc(flipud(rs3(1:128,:)))
set(gca,'YTick',[]);set(gca,'XTick',[])
xlabel('time')
ylabel('frequency')
title('mode #3')