關(guān)于Hilbert-Huang的matlab實(shí)現(xiàn)，材料匯總，比較雜...感謝所有網(wǎng)絡(luò)上的貢獻(xiàn)者們:)
核心：以下代碼計(jì)算HHT邊際譜及其對(duì)應(yīng)頻率
工具包要求：G-Rilling EMD Toolbox，TFTB Toolbox
附：黃鍔先生課題組開(kāi)發(fā)的工具包（可以在 這里 找到），這里并未用到。
% Empirical mode decomposition, resulting in intrinc mode functions.% Without parameter 'MAXMODES' the process may be seriously delayed by% decompose original signals into too many IMFs (not necessary, 9 is % enough generally)imfs = emd(oriSig, 'MAXMODES', 9);% HHT spectrum: hhtS[A, f, t] = hhspectrum(imfs);[hhtS, ~, fCent] = toimage(A, f, t);% Marginal hilbert spectrum: hhtMS, xf: correspondig frequencyfor k = 1 : size(hhtS, 1) hhtMS(k) = sum(hhtS(k, : )) * 1 / fs;endxf = fCent(1, :) .* fs;簡(jiǎn)單來(lái)說(shuō)，設(shè)置好路徑之后輸入 install_emd 即可。
% [A,f,tt] = HHSPECTRUM(imf,t,l,aff)% Input:%- imf : matrix with one IMF per row % 將emd分解得到的IMF代入就可以，就是你的程序中寫(xiě)的c變量，不用加最后一行的趨勢(shì)項(xiàng)%- t : time instants % 瞬時(shí)時(shí)間或持續(xù)時(shí)間 ？？（寫(xiě)[1:信號(hào)長(zhǎng)度]就可以，真實(shí)的時(shí)間可以根據(jù)采樣率轉(zhuǎn)換%- l : estimation parameter for instfreq % 瞬時(shí)頻率的估計(jì)參數(shù) ？？（寫(xiě)1就可以，決定瞬時(shí)頻率估計(jì)時(shí)的邊界從第幾個(gè)點(diǎn)開(kāi)始%- aff : if 1, displays the computation evolution % 顯示計(jì)算進(jìn)程選項(xiàng)，不想顯示寫(xiě)0就可以%% Output:% - A : amplitudes of IMF rows% - f : instantaneous frequencies% - tt : truncated time instants % 截止時(shí)間 ？？（截?cái)鄷r(shí)間，返回的是瞬時(shí)頻率對(duì)應(yīng)的時(shí)間，要比原來(lái)信號(hào)的時(shí)間按短，由前面的l值決定）%% calls:% - hilbert : computes the analytic signal% - instfreq : computes the instantaneous frequency % 瞬時(shí)頻率%% Example：[A,f,tt] = hhspectrum(c(1:end-1,:),[1:n],1,0);[im,tt,ff] = toimage(A,f,tt,512);disp_hhs(im,tt,[],fs);關(guān)于時(shí)頻圖的概念，我認(rèn)為是與諸如小波，gabor等聯(lián)合時(shí)頻分析方法聯(lián)系在一起的。
小波，gabor等具有多尺度分析的概念，得到的時(shí)頻分布是一個(gè)二維的矩陣（橫軸時(shí)間，縱軸頻率，可以用不同的顏色（光譜圖spectrogram）或瀑布圖表示不同的幅度）。
對(duì)HHT來(lái)講，我覺(jué)得氣時(shí)頻圖的概念是有些不同的。EMD分解的作用就是把復(fù)雜的信號(hào)分界為簡(jiǎn)單的單分量的信號(hào)，使其可以應(yīng)用瞬時(shí)頻率的概念，hilbert變換的目的就是分析出瞬時(shí)頻率。所以HHT在每一時(shí)刻得到的只有一個(gè)值，而不是像小波之類(lèi)的得到一系列的值（多尺度分析）。所以我們從其時(shí)頻分布圖上看到的是一條線，而不是一幅圖。
%示例程序N=1000;n=1:N;fs=1000;t=n/fs;fx=50;fy=150;x=2*sin(2*pi*fx*t);y=sin(2*pi*fy*t);z=x+y;data=z;imf=emd(data); % 對(duì)輸入信號(hào)進(jìn)行EMD分解 [A,f,t]=hhspectrum(imf); % 對(duì)IMF分量求取瞬時(shí)頻率與振幅：A：是每個(gè)IMF的振幅向量,f:每個(gè)IMF對(duì)應(yīng)的瞬時(shí)頻率，t：時(shí)間序列號(hào)[E,t,Cenf]=toimage(A,f); % 將每個(gè)IMF信號(hào)合成求取Hilbert譜，E：對(duì)應(yīng)的振幅值，Cenf：每個(gè)網(wǎng)格對(duì)應(yīng)的中心頻率。這里橫軸為時(shí)間，縱軸為頻率。即時(shí)頻圖(用顏色表示第三維值的大小)和三維圖（三維坐標(biāo)系：時(shí)間，中心頻率，振幅） cemd_visu(data,1:length(data),imf); % 顯示每個(gè)IMF分量及殘余信號(hào)disp_hhs(E); % 希爾伯特譜% 畫(huà)出邊際譜colormap(flipud(gray)); % 黑白顯示%N=length(Cenf); % 設(shè)置頻率點(diǎn)數(shù)。完全從理論公式出發(fā)。網(wǎng)格化后中心頻率很重要，大家從連續(xù)數(shù)據(jù)變?yōu)殡x散的角度去思考，相信應(yīng)該很容易理解for k=1:size(E,1) bjp(k)=sum(E(k,:))*1/fs;endfigure(3);plot(Cenf(1,:)*fs,bjp); % 作邊際譜圖 進(jìn)行求取Hilbert譜時(shí)頻率已經(jīng)被抽樣成具有一定窗長(zhǎng)的離散頻率，所以此時(shí)的頻率軸已經(jīng)是中心頻率xlabel('頻率 / Hz');ylabel('幅值');% 繪制瞬時(shí)包絡(luò)圖和瞬時(shí)頻率圖figure;subplot(221),plot(t/N,A(1,:));xlabel('時(shí)間 t/s');ylabel('幅值');title('imf1分量瞬時(shí)包絡(luò)');subplot(222),plot(t/N,f(1,:)*fs);xlabel('時(shí)間 t/s');ylabel('頻率');title('imf1分量瞬時(shí)頻率');subplot(223),plot(t/N,A(2,:));xlabel('時(shí)間 t/s');ylabel('幅值');title('imf2分量瞬時(shí)包絡(luò)');subplot(224),plot(t/N,f(2,:)*fs);xlabel('時(shí)間 t/s');ylabel('頻率');title('imf2分量瞬時(shí)頻率');Hilbert譜：信號(hào)的希爾伯特變換后做fft，表示信號(hào)幅值在整個(gè)頻率段上隨時(shí)間和頻率的變化規(guī)律；
Hilbert邊際譜：對(duì)hilbert譜做積分，表示信號(hào)幅值在整個(gè)頻率段上隨頻率的變化情況，它相當(dāng)于傅里葉譜，但比傅里葉譜具有更高的頻率分辨率。Hilbert邊際譜是通過(guò)對(duì)Hilbert譜（在時(shí)間上）積分得到的；
Hilbert包絡(luò)譜：希爾伯特變換后做包絡(luò)后再fft，不同于Hilbert譜和Hilbert邊際譜，是直接對(duì)信號(hào)進(jìn)行Hilbert變換后構(gòu)造解析函數(shù)，然后依據(jù)解析函數(shù)求模值，求的模值即為包絡(luò)，然后對(duì)信號(hào)包絡(luò)進(jìn)行FFT后得到的即為Hilbert包絡(luò)譜。
邊際譜從統(tǒng)計(jì)意義上表示了整組數(shù)據(jù)每個(gè)頻率點(diǎn)的積累幅值分布，而傅立葉譜的某點(diǎn)幅值表示在整個(gè)信號(hào)里有一個(gè)含有此頻率的三角函數(shù)組分，而且幅值越大只是說(shuō)明在整個(gè)數(shù)據(jù)段上，局部存在的可能性越大。
再看得到的圖形，F(xiàn)FT 表示的是整個(gè)數(shù)據(jù)中，能量在一個(gè)頻率上分布的可能性地描述，而邊際譜表示在在每一個(gè)頻率上幅值的積累，如果想知道具體時(shí)間那么就看HHT譜，這個(gè)時(shí)間-幅值-頻率的三維譜。說(shuō)到瞬時(shí)頻率，傅立葉變換不強(qiáng)調(diào)局部性，而是強(qiáng)調(diào)全局性。咱們的HHT才提出一個(gè)唯一的瞬時(shí)頻率的定義。因此拿瞬時(shí)頻率來(lái)衡量傅立葉變換也是不公平的。
對(duì)于邊際譜，就是hilbert譜對(duì)時(shí)間的積分，從積分的角度來(lái)講，就相當(dāng)于對(duì)任意一階頻率把所有的時(shí)間上的幅值都加起來(lái)了，這就反映這階頻率在所有時(shí)間上的幅值積累，而對(duì)于他們對(duì)于頻率所對(duì)應(yīng)的幅值，描述為該頻率在整個(gè)時(shí)間上出現(xiàn)的可能，我個(gè)人認(rèn)為是既然出現(xiàn)了，那就是存在的，不能說(shuō)是一種可能，只能說(shuō)他出現(xiàn)的次數(shù)的多少，經(jīng)過(guò)累積以后就變成了他的幅值，而對(duì)于fourier來(lái)講，他只能說(shuō)明在用正余弦函數(shù)擬合這個(gè)信號(hào)的時(shí)候需要這一階頻率，幅值大，則說(shuō)明他對(duì)擬合這個(gè)信號(hào)的貢獻(xiàn)大，而不是一種出現(xiàn)的可能性的度量。
能量譜，可以理解為對(duì)邊際譜的平方，這個(gè)只是具有能量的形式，而具體是否能代表能量，需要更進(jìn)一步的探討，還有一個(gè)瞬時(shí)能量，這些名詞確實(shí)很誘人，但究竟如何，還需要大家的努力。
聯(lián)系客服
