% EMD程序中文詳解及應(yīng)用例子
% 
% 
% 看到版上總是有人問(wèn)EMD程序應(yīng)用方面的問(wèn)題，于是把 G.Rilling 寫的 EMD.m注釋漢化了，同時(shí)根據(jù)個(gè)人理解添加了部分注釋，最后附上一個(gè)EMD分解＋HHT譜＋邊際譜的例子。例子程序中所用到的函數(shù)是G.Rilling提供的，可以到如下地址下載：
% http://perso.ens-lyon.fr/patrick.flandrin/emd.html 
% http://zhlong.ys168.com/
% 由于個(gè)人水平有限，權(quán)作拋磚引玉，希望對(duì)各位研究EMD的朋友有所幫助，也歡迎大家指正。

% 計(jì)算2FSK信號(hào)的HHT譜和邊際譜
% 作者：xray    2007.11
clear

rand('seed', 0);

T = 0.05;                % 仿真時(shí)間
R = 500;                 % 碼速率
fd = 1000;               % 載波頻差
fc = 2000;               % 載波頻率
fs = 20000;              % 采樣率
samp = fs/R;             % 每個(gè)碼元上的采樣點(diǎn)數(shù)
N = T*fs;
n = 1:N;
x = randint(1, R*T, 2);
y = fskmod(x, 2, fd, samp, fs);
y = y .* exp(i*2*pi*fc/fs*n);
y = real(y);
% z = awgn(y, 20, 'measured'); 
z = y;

imf = emd(z); 
[A, fa, tt] = hhspectrum(imf);
if size(imf,1) > 1
    [A,fa,tt] = hhspectrum(imf(1:end-1, :));
else
    [A,fa,tt] = hhspectrum(imf);
end
[E, tt1] = toimage(A,fa,tt,length(tt));    

for k = 1:size(E,1)
    bjp(k) = sum(E(k,:))*1/fs*1/T;  
end
f = (0:N-3)/N*(fs/2);

figure(1)
plot(z);

figure(2)
imagesc(tt1,[0,0.5*fs],E);
set(gca,'YDir','normal')
% 使用灰度圖顯示
% colormap(flipud(gray))

figure(3)
plot(f, bjp);

% 例子說(shuō)明：
% (1) 運(yùn)行該實(shí)例還需要時(shí)頻工具箱(tftb工具包)支持，可以到　http://zhlong.ys168.com/　下載
% (2) 由于該程序調(diào)用的m文件位于多個(gè)文件夾，因此最好將如下路徑
% *\emd\emds
% *\emd\utils
% *\tftb-0.1\mfiles
% 添加到matlab的路徑中(File\Set Path)
% (3) 2FSK是通信中常用的數(shù)字調(diào)制方式，分別用兩個(gè)頻率的正弦波表示“1”和“0”，例子中這兩個(gè)頻率分別為1500Hz和2500Hz。圖1是2FSK信號(hào)的時(shí)域波形，可以看到振動(dòng)波形有疏有密。對(duì)照?qǐng)D2可以發(fā)現(xiàn)，疏密部分對(duì)應(yīng)了信號(hào)中的低頻部分(1500Hz)和高頻部分(2500Hz)，HHT譜恰好實(shí)現(xiàn)了2FSK信號(hào)的解調(diào)，這也是這個(gè)應(yīng)用的物理意義。圖3反映了信號(hào)在兩個(gè)頻率上的能量分布，表示了信號(hào)中“1”和“0”中出現(xiàn)的比例。