1. 关于经验模态分解(EMD)的问题
于非线性时间序列分析经验模态分解和小波分解异同性的研究
龚志强 邹明玮 高新全 董文杰
摘 要：基于经验模态分解(EMD)的希尔伯特变换(HT),是对非线性时间序列基于EMD进行分解,然后通过HT获得频谱.利用理想时间序列和青藏高原古里雅冰芯18O时间序列,系统地分析比较了EMD和小波分解(WD)以及HT和小波变换在非线性时间序列处理中的优劣,并针对它们各自的缺点提出了可能改进的设想.研究结果表明,将基于EMD的方法和基于WD的方法有机结合起来应用,可以更有效地识别原时间序列的特征信息.
关键词：经验模态分解;小波分解;理想时间序列;古里雅冰芯
文章编号：1000-3290/2005/54(08)/3947-11
On the difference between empirical mode decomposition and wavelet decomposition in the nonlinear time series
Gong Zhi-Qiang Zou Ming-Wei Gao Xin-Quan Dong Wen-Jie
基金项目：国家重点基础研究发展规划(批准号:2004CB418300)和国家自然科学基金(批准号:90411008,40231006)资助的课题.
作者单位：龚志强（扬州大学物理科学与技术学院,扬州,225009;国家气候中心气候研究开放实验室,北京,100081）
邹明玮（扬州大学物理科学与技术学院,扬州,225009;中国科学院大气物理研究所,北京,100029）
高新全（国家气候中心气候研究开放实验室,北京,100081）
董文杰（国家气候中心气候研究开放实验室,北京,100081）
参考文献：
[1]Huang N E, Shen Z, Long R S et al1998 Proc. R. Soc. Lond. A454 899
[2]Huang N E 1999 Ann. Rev. Fluid Mech. 31 417
[3]Feng G L, Chou J M, Dong W J 2004 Chin. Phys. 13 1582
[4]Feng G L, Dong W J 2003 Chin. Phys. 13 413
[5]Dai X G, Wang P, Chou J F 2004 Prog. Nat. Sci. 14 73
[6]Li J P, Tang Y Y 1999 Application of Wavelet Analysis ( Chongqing:Chongqing University Press)(in Chinese)[李建平、唐远延1999小波分析方法的应用(重庆:重庆大学出版社)]
[7]Mallat S 1989 IEEE Trans. Signal Proces. 37 2091
[8]Mallat S, Hwang W L 1992IEEE Trans. Inform. Theory 38 617
[9]Xiong X J, Guo B H, Xu Y M 2002 J. Ocean. HB Seas 20 12(in Chinese)[熊学军、郭炳火、胡筱敏2002黄渤海海洋20 12]
[10]Farge M 1992 Ann. Rev. Fluid Mech. 24 395
[11]Tewfiki A H 1992 IEEE Trans. Inform. Theory 38 747
[12]Li B B 1994 J. Electron. 16 646(in Chinese)[李兵兵 1994 电子科学学刊16 646]
[13]Wang J Z 1998 Auto.Elec.Power Sys.22 40(in Chinese)[王建赜1998电力系统自动化22 40]
[14]Yao T D,Yang X M,Kang X C 2001 Quatern.Sci.21 514(in Chinese)[姚檀栋、杨学梅、康兴成2001第四纪研究21 514]
[15]Zhang X P,Yao T D,Jin H J 2000 J.Glaci.Geocr.22 23(in Chinese)[章新平、姚檀栋、金会军2000冰川冻土22 23]
[16]Zhang X P,Shi Y F,Yao T D 1995 Sci.China D 38 854(in Chinese)[章新平、施雅风、姚檀栋1995中国科学D 38 854]
[17]Yao T D,Jiao K Q,Li Z Q 1994 Sci.China B 24 763(in Chinese)[姚檀栋、焦克勤、李忠勤1994中国科学B 24 763]
[18]Briffa K R 1998 Nature 391 678
[19]Briffa K R 1998 Nature 393 450
[20]Stagke D W 1998 Science 280 564
[21]Briffa K R 1998 Quatern. Sci. Rev. 19 87
[22]Liu C X2001 J.Trop.Meteor.17 381(in Chinese)[刘春霞2001热带气象学报17 381]
2. linux环境c语言编程书籍，那种从初级入门到进阶深层次一点的
其实APUE这本书很好，灰常之经典，不知道啥意思？就是LS那位仁兄说的《unix环境高级编程》，英文名是Advanced Programming In Unix Environment ， 建议你还是阅读英文版，计算机类图书生词都很少，查查字典就能看懂了，中文版翻译的不给力，还有就是不要觉得这本书难，只有起点高进步才会快，看不懂就反复的读。
3. 可以分享一下基于emd的小波去噪程序吗谢谢啦。
供参考：
lev=5;
[c,l]=wavedec(x,lev,wname);
sigma=wnoisest(c,l,1);
alpha=2;
thr1=wbmpen(c,l,sigma,alpha)
[thr2,nkeep]=wdcbm(c,l,alpha)
xd1=wdencmp('gbl',c,l,wname,lev,thr1,'s',1);
[xd2,cxd,lxd,perf0,perfl2]=wdencmp('lvd',c,l,wname,lev,thr2,'h');
[thr,sorh,keepapp]=ddencmp('den','wv',x)
xd3=wdencmp('gbl',c,l,wname,lev,thr,'s',1);
subplot(411);plot(x);title('原始信号','fontsize',12);
subplot(412);plot(xd1);title('使用penalty阈值降噪后信号','fontsize',12);
subplot(413);plot(xd2);title('使用Birge-Massart阈值降噪后信号','fontsize',12);
subplot(414);plot(xd3);title('使用缺省阈值降噪后信号','fontsize',12);
s=[-1.58 0.42 0.46 0.78 -0.49 0.59 -1.3 -1.42 -0.16 -1.47 -1.35 0.36 -0.44 -0.14 1 -0.5 -0.2 -0.06 -0.6 0.42 -1.52 0.51 0.76 -1.5 0.16 -1.29 -0.65 -1.48 0.6 -1.65 -0.55];
[C,L]=wavedec(s,1,'db3');
ca1=wrcoef('a',C,L,'db3',1);
x1=ca1 ;
[C,L]=wavedec(s,2,'db3');
ca2=wrcoef('a',C,L,'db3',2);
x2=ca2 ;
[C,L]=wavedec(s,3,'db3');
ca3=wrcoef('a',C,L,'db3',3);
x3=ca3 ;
[C,L]=wavedec(s,4,'db3');
ca4=wrcoef('a',C,L,'db3',4);
x4=ca4 ;
cg = wrcoef('a',C,L,'sym5',1);
x5=cg;
p=1:31;
subplot(6,1,1);plot(p,s);ylabel('s');
subplot(6,1,2);plot(p,x1);ylabel('ca1');
subplot(6,1,3);plot(p,x2);ylabel('ca2');
subplot(6,1,4);plot(p,x3);ylabel('ca3');
subplot(6,1,5);plot(p,x4);ylabel('ca4')
subplot(6,1,6);plot(p,x5);ylabel('ca5') %加入的重构，是不是你要的？
4. 求助，请教EMD工具箱使用问题，出来的c2f，f2c是什么东西
你在命令床后输入“help emd”就知道这个工具箱是否添加进来了，没有的话就是路径设置不成功，可能是没有保存
5. 求大神给我帮个忙。matlab的s-function怎么调用不了emd（文件）函数啊看下面的补充问题。。万分感谢。
你的问题可以化为下面向量的问题
已知a=(1,1,1)，b=(-1,1,1)，c=a×u，d=c×u，
c和d的夹角是50°，c和v的夹角是55°，d和v的夹角是4.9°，
u⊥v，|u|=1，|v|=1
求u，v
题中的a，b,c,d,u,v均为三维向量，×表示向量内积，|u|表示向量u的模
其中，向量b对应你以前的（m,n,p），向量u对应你以前的（h,k,l）,向量v对应你以前的(u,v,w)
由上题，c=a×u，d=c×u可得c⊥u，d⊥u又u⊥v，且c,d,v有相同的起点即坐标原点，从而c,d,v在同一平面上且有相同的起点，且均与u垂直
所以c,d,v之间的夹角必定满足某个等式，回到题上也就是说，55°=50°+5°，
进一步说，你给的条件是矛盾的，所以matlab找不到解
就算你给出的条件是对的，由于你给出的前三个方程并非完全独立的，也不足以确定你想要的结果
追问
这是其中的一部分条件，我就是想问一下，怎么去编写一个不确定的参数，让它可以在矩阵中使用，还不出现Inner matrix dimensions must agree.等错误
6. 谁可以给我一个emd分解的matlab程序。只需要emd分解的。谢谢了！
%此版本为ALAN 版本的整合注释版
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为极值点所在位置
% 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:10; % 采样率2000HZ
% 调幅信号
%x=sin(2*pi*t).*sin(40*pi*t);
x=sin(2*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
7. 在运行里面 输入“emd”后，出现的那个对话框是什么东西
应该是cmd
cmd是command的缩写.命令行
在9x系统下输入command就可以打开命令行.而在NT系统上可以输入cmd来打开.
操作顺序是：开始－＞运行－＞键入cmd或command
在命令行里你可以看到你的系统版本,文件系统版本等等
你可以敲入help查看帮助
Cmd启动命令解释器Cmd.exe的新实例。如果在不含参数的情况下使用，则cmd显示WindowsXP的版本和版权信息。
语法
cmd[[{/c|/k}][/s][/q][/d][{/a|/u}][/t:fg][/e:{on|off}][/f:{on|off}][/v:{on|off}]string]
参数
/c
执行string指定的命令，然后停止。
/k
执行string指定的命令并继续。
/s
修改位于/c或/k之后的string处理。
/q
关闭回显。
/d
禁用自动运行命令执行。
/a
创建美国国家标准协会(ANSI)输出。
/u
创建Unicode输出。
/t:fg
设置前景f和背景g的颜色。下表列出了可用作f和g的值的有效十六进制数字。值颜色
0黑色
1蓝色
2绿
3湖蓝色
4红
5紫色
6黄
7白色
8灰色
9浅蓝色
A浅绿色
B浅水绿
C浅红色
D浅紫色
E浅黄色
F亮白色
/e:on
启用命令扩展。
/e:off
禁用命令扩展。
/f:on
启用文件和目录名完成。
/f:off
禁用文件和目录名完成。
/v:on
启用延迟的环境变量扩展。
/v:off
禁用延迟的环境变量扩展。
string
指定要执行的命令。
/?
在命令提示符显示帮助。
注释
使用多个命令
可以在string中使用由&&分隔的多个命令，不过这些命令必须置于引号之中（例如，"command&&command&&command"）。
处理引号
如果指定了/c或/k，则在满足下述所有条件的情况下，cmd会处理string中的其余命令而将引号保留：
未使用/s。
正确使用一对引号。
在引号内未使用任何特殊字符（例如：&<>()@^|}。
在引号内使用了一个或多个空格子符。
引号内的string为可执行文件的名称。
如果上述条件不能满足，则处理string时将首先检查它的第一个字符以验证其是否为左引号。如果第一个字符是左引号，则它会与右引号分离开。跟在右引号之后的任何文本都会得到保留。
执行注册表子项
如果在string中未指定/d，Cmd.exe会查找下述注册表子项：
HKEY_LOCAL_MACHINE\Software\Microsoft\CommandProcessor\AutoRun\REG_SZ
HKEY_CURRENT_USER\Software\Microsoft\CommandProcessor\AutoRunREG_EXPAND_SZ
如果上述的一个注册表子项或两个都存在，则会在执行其他变量之前执行它们。
警告
编辑注册表不当可能会严重损坏您的系统。在更改注册表之前，应备份计算机上任何有价值的数据。
启用和禁用命令扩展
在WindowsXP中，命令扩展在默认情况下是启用的。对于特定过程可以使用/e:off将它们禁用。通过设置下述REG_DWORD值，可以在计算机上或用户会话中启用或禁用所有cmd命令行选项的扩展：
HKEY_LOCAL_MACHINE\Software\Microsoft\CommandProcessor\EnableExtensions\REG_DWORD
HKEY_CURRENT_USER\Software\Microsoft\CommandProcessor\EnableExtensions\REG_DWORD
在注册表中使用Regedit.exe可以将REG_DWORD值设为0×1（即启用）或0×0（即禁用）。用户特定设置优先于计算机设置，并且命令行选项优先于注册表设置。
警告
编辑注册表不当可能会严重损坏您的系统。在更改注册表之前，应备份计算机上任何有价值的数据。
启用命令扩展后，会影响到下述命令：
assoc
call
chdir(cd)
color
del(erase)
endlocal
for
ftype
goto
if
mkdir(md)
popd
prompt
pushd
set
setlocal
shift
start（还包括将更改外部命令过程）
有关这些命令的详细信息，请参阅“相关主题”。
启用延迟的环境变量扩展
启用延迟的环境变量扩展，可以使用感叹号字符来替代运行时的环境变量值。
启用文件和目录名完成
默认情况下，禁用文件和目录名完成。对于特定的cmd命令处理，可以通过/f:{on|off}来启用或禁用该功能。通过设置下述REG_DWORD值，可以在计算机上或用户会话中启用或禁用所有cmd命令处理的文件和目录名完成：
HKEY_LOCAL_MACHINE\Software\Microsoft\CommandProcessor\CompletionChar\REG_DWORD
HKEY_LOCAL_MACHINE\Software\Microsoft\CommandProcessor\PathCompletionChar\REG_DWORD
HKEY_CURRENT_USER\Software\Microsoft\CommandProcessor\CompletionChar\REG_DWORD
HKEY_CURRENT_USER\Software\Microsoft\CommandProcessor\PathCompletionChar\REG_DWORD
要设置REG_DWORD值，请运行Regedit.exe并使用特定功能的控制字符的十六进制值（例如，用0×9表示TAB键，用0×08表示BACKSPACE键）。用户特定设置优先于计算机设置，并且命令行选项优先于注册表设置。
警告
编辑注册表不当可能会严重损坏您的系统。在更改注册表之前，应备份计算机上任何有价值的数据。
如果使用/f:on启用了文件和目录名完成，则对于目录名完成，可使用CTRL+D组合键；而对于文件名完成，可使用CTRL+F组合键。要禁用注册表中特定字符的完成，请使用空格值[0×20]，因为空格不是有效的控制字符。
按CTRL+D或CTRL+F组合键时，cmd会处理文件和目录名的完成操作。这些组合键的作用是在string后附加通配符（如果还未使用），并创建匹配的路径列表，然后显示第一个匹配的路径。如果所有路径都不匹配，文件和目录名完成操作会发出警告声，并且不更改所显示的内容。要逐个查看匹配路径列表中的路径，请重复按CTRL+D或CTRL+F组合键。要向后查看该列表，请在按SHIFT的同时按CTRL+D或CTRL+F组合键。要放弃已保存的匹配路径列表并生成新列表，可以编辑string，然后按CTRL+D或CTRL+F组合键。如果在CTRL+D和CTRL+F组合键之间切换，将会放弃已保存的匹配路径列表并生成新列表。CTRL+D组合键与CTRL+F组合键之间唯一的不同在于，CTRL+D仅匹配目录名，而CTRL+F既匹配文件名，又匹配目录名。如果在任何内部目录命令（CD、MD或RD）中使用文件和目录名的完成，将仅使用目录的完成。
如果将匹配路径置于引号之中，则文件和目录名完成会正确地处理含有空格或特殊字符的文件名。
下述特殊字符需要有引号：&<>[]{}^=;!'+,`~[whitespace]
如果您提供的信息包含空格，请将文本置于引号之中（例如，"ComputerName"）。
如果从string中处理文件和目录名完成操作，则位于光标右侧的[Path]的任意部分都将放弃（即在string中处理完成操作的位置）。
格式化图例
格式含义
斜体用户必须提供的信息
粗体用户必须像显示的一样准确键入的元素
省略号(...)可在命令行中重复多次的参数
在括号([])之间可选项目
在大括号({})之间；将选项用竖线(|)隔开。例如：{even|odd}用户必须从中只选择一个选项的选项组
Courier字体代码或程序输出
8. matlab装上EMD工具箱后怎么调用emd函数啊
下载的工具箱中的emd.m文件里的注释有详细的用法介绍。
工具箱的安装
运行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环境中运行，算法的主要功能函数并不受影响。
9. 怎么用c语言写状态机呀请举例说明
c语言写状态机之前：
1、确定一共有多少种状态,这里的状态有开和关,细分还有say thankyou 和警告
2、确定状态之间的迁移条件
如果按照四种状态：开、关、说谢谢、警告,那么这四种状态之前的迁移条件很明显了
分两个函数：
1、检查是否需要迁移状态；
2、迁移状态.
遍历各种状态检查是否有状态需要发生迁移.一般用一个switch将各种状态列出,然后在各种状态里面用if检查是否需要迁移状态,如果需要迁移,做好标记.
再次遍历各种状态,检查哪些状态做了标记,迁移到新状态,并做相应的操作,比如进入关的时候,做关门动作。
典型的状态机结构：
enum { state_A, state_B, state_C } state = state_A;
while(1)
{
switch ( state )
{
case state_A:
if ( event_A ) // 这里也可以用switch
{
action_1(); // 在某状态下发生某事件执行某个动作，并转入下个状态
state = state_B;
}
else if ( event_B )
{
}
else
{
}
break;
case state_B:
... ...
}
}
10. 跪求一个emd 去噪的程序 matlab 代码 带中文解释的 方便理解
function imf = emd(x,n);%%最好把函数名改为emd1之类的，以免和Grilling的emd冲突
%%n为你想得到的IMF的个数
c = x('; % of the input signal (as a row vector)
N = length(x);-
% loop to decompose the input signal into n successive IMFs
imf = []; % Matrix which will contain the successive IMF, and the resiefor t=1:n
% loop on successive IMFs
%-------------------------------------------------------------------------
% inner loop to find each imf
h = c; % at the beginning of the sifting process, h is the signal
SD = 1; % Standard deviation which will be used to stop the sifting process
while SD > 0.3 % while the standard deviation is higher than 0.3 (typical value) %%筛选停止准则
% find local max/min points
d = diff(h); % approximate derivative %%求各点导数
maxmin = []; % to store the optima (min and max without distinction so far)
for i=1:N-2
if d(i)==0 % we are on a zero %%导数为0的点，即”驻点“，但驻点不一定都是极值点，如y=x^3的x=0处
if sign(d(i-1))~=sign(d(i+1)) % it is a maximum %%如果驻点两侧的导数异号（如一边正，一边负），那么该点为极值点
maxmin = [maxmin, i]; %%找到极值点在信号中的坐标（不分极大值和极小值点）
end
elseif sign(d(i))~=sign(d(i+1)) % we are straddling a zero so%%如y=|x|在x=0处是极值点，但该点倒数不存在，所以不能用上面的判
断方法
maxmin = [maxmin, i+1]; % define zero as at i+1 (not i) %%这里提供了另一类极值点的判断方法
end
end
if size(maxmin,2) < 2 % then it is the resie %%判断信号是不是已经符合残余分量定义
break
end
% divide maxmin into maxes and mins %% 分离极大值点和极小值点
if maxmin(1)>maxmin(2) % first one is a max not a min
maxes = maxmin(1:2:length(maxmin));
mins = maxmin(2:2:length(maxmin));
else % is the other way around
maxes = maxmin(2:2:length(maxmin));
mins = maxmin(1:2:length(maxmin));
end % make endpoints both maxes and mins
maxes = [1 maxes N];
mins = [1 mins N];
%------------------------------------------------------------------------- % spline interpolate to get max and min envelopes; form imf
maxenv = spline(maxes,h(maxes),1:N); %%用样条函数插值拟合所有的极大值点
minenv = spline(mins, h(mins),1:N); %%用样条函数插值拟合所有的极小值点
m = (maxenv + minenv)/2; % mean of max and min enveloppes %%求上下包络的均值
prevh = h; % of the previous value of h before modifying it %%h为分解前的信号
h = h - m; % substract mean to h %% 减去包络均值
% calculate standard deviation
eps = 0.0000001; % to avoid zero values
SD = sum ( ((prevh - h).^2) ./ (prevh.^2 + eps) ); %% 计算停止准则
end
imf = [imf; h]; % store the extracted IMF in the matrix imf
% if size(maxmin,2)<2, then h is the resie
% stop criterion of the algo. if we reach the end before n
if size(maxmin,2) < 2
break
end
c = c - h; % substract the extracted IMF from the signal
end
return
参考资料
http://..com/link?url=5jhgpTXu95SDRgi_