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_