MATLAB實作EMD分解及希爾伯特譜分析

希爾伯特—黃變換

傳統的傅里葉變化只能得到信號在采樣周期內的全域頻率資訊，處理頻率隨時間變化的非平穩信號具有很大的局限性，希爾伯特－黃變換(Hilbert-Huang Trans form) 是由N. E. Huang 等人于1998年提出的一種非線性、非平穩信號的分析處理方法，可以求得信號的瞬時頻率、瞬時幅值等瞬時特征，

原理

希爾伯特變換：

g ∧ ( t ) = 1 π P ∫ ? ∞ + ∞ g ( τ ) t ? τ d τ = P π t ? g ( t ) \mathop g\limits^ \wedge (t) = \frac{1}{\pi }P\int\limits_{ - \infty }^{ + \infty } {\frac{{g(\tau )}}{{t - \tau }}} d\tau = \frac{P}{{\pi t}} * g(t) g∧?(t)=π1?P?∞∫+∞?t?τg(τ)?dτ=πtP??g(t)

信號g(t) 的希爾伯特變換是原信號g(t) 與1/ π t \pi t πt) 在時域內的卷積，這個卷積在時域內很難理解；將它轉換到頻域中來看，時域中的卷積相當于頻域中的相乘，因此，希爾伯特變換的結果是給原來的實信號提供一個幅值和頻率不變，但相位平移90° 的信號，希爾伯特變換將原信號的延遲90°相位信號作為虛部，原信號作為實部，構成一個復數信號，從而在復數域求解信號的瞬時特征，如下圖：

利用matlab自帶的hht函式畫“時間—頻率—幅值圖”(無背景色）

fs=500;
f1=15;
tmax=1;
N=tmax*fs;     %采樣點數
t=0:1/fs:tmax-1/fs;
x=sin(t*2*pi*f1);
hht(x,fs,'FrequencyLimits',[0,50])   %畫二維顏色圖

自創函式畫有背景色的“時間—頻率—幅值圖”

fs=500;
f1=15;
tmax=1;
N=tmax*fs;     %采樣點數
t=0:1/fs:tmax-1/fs;
x=sin(t*2*pi*f1);
H_x=hilbert(x);         %希爾伯特變換
df=abs(diff(atan(imag(H_x)./x0))./(2*pi*diff(t)));   %差分求瞬時頻率
E=abs(H_x);             %求瞬時幅值
tt=t(1:length(df));
df_c=fs/(10*N);
dtt_c=tmax/N;
%轉換為可以畫RGB圖的矩陣
for i=1:N
    for j=1:N-1
        if((df_c*i)>df(j)&&(df_c*(i-1))<df(j))
            C(i,j)=E(j);
        end
    end
end
imagesc([0 tmax],[0 fs/10],C)         %畫RGB圖
xlabel('t'); 
ylabel('f');
title('hht')
colorbar

這里由于imagesc函式y軸的值是從上往下遞增的，所以影像和hht畫的圖相反

EMD（經驗模態分解）

瞬時頻率對信號的約束條件，給人們一種啟示：在對信號進行希爾伯特變換之前，要先把信號分解為瞬時頻率具有意義的各個分量，把信號進行分解的方法在希爾伯特－黃變換理論中稱為經驗模式分解，在此方法中N. E. Huang 等人定義了一類函式，稱為固有模式函式，利用這類函式的區域特性，可以使函式在任意一點的瞬時頻率都有意義，經驗模式分解的最大貢獻是使信號符合了“單分量”的要求，進而使瞬時頻率有了明確的物理意義，從而使得希爾伯特－黃變換比單純的希爾伯特變換在時頻分析方面有了很大地進步，

具體原理就不詳述了，直接給出代碼~

EMD分解后各階IMF（固有模態分量）的“時間—頻率—幅值圖”

fs=500;
f1=15;
f2=10;
tmax=1;
N=tmax*fs;     %采樣點數
t=0:1/fs:tmax-1/fs;
noise=wgn(1,N,2);
x=sin(t*2*pi*f1);%+sin(t*2*pi*f2);
y=x+noise;

figure(1);
subplot(211);
plot(t,y);
title('Original Signal');
xlabel('Time'); 
ylabel('Amplitude'); 
subplot(212);
[fft_f,fft_a]=fftlw(t,y);
plot(fft_f,fft_a);
title('Original Signal');
xlabel('F'); 
ylabel('Amplitude');

figure(2)
[imf,r]=emd(y);
for i=1:3
    subplot(410+i)
    plot(t,imf(:,i))
end
subplot(414)
plot(t,r)


x0=imf(:,3)';
H_x=hilbert(x0);
figure(3);
subplot(311);
plot(t,x0);
title('Original imf');
xlabel('t'); 
ylabel('imf'); 
subplot(312);
plot(t,imag(H_x));
xlabel('t'); 
ylabel('imag(imf)');
subplot(313);
plot(x0,imag(H_x));
xlabel('imf'); 
ylabel('imag(H_imf)');

C=zeros(N,N);
for p=1:length(imf(1,:))
    x0=imf(:,p)';
    H_x=hilbert(x0);
    df=abs(diff(atan(imag(H_x)./x0))./(2*pi*diff(t)));
    %df=gradient(atan(imag(H_x')./x0))./(gradient(t)*2*pi);
    E=abs(H_x);
    tt=t(1:length(df));
    df_c=fs/(10*N);
    dtt_c=tmax/N;
    for i=1:N
        for j=1:N-1
            if((df_c*i)>df(j)&&(df_c*(i-1))<df(j))
                C(i,j)=E(j);
            end
        end
    end
end


            

figure(4)
imagesc([0 tmax],[0 fs/10],C)
xlabel('t'); 
ylabel('f');
title('hht')
% imagesc(tt,df,E)
colorbar

figure(5)
hht(imf,fs,'FrequencyLimits',[0 50])

經對比可以發現，自創函式做出的圖（需要倒過來一下）與matlab自帶的hht函式相比，matlab自帶的hht函式做出的圖較雜亂~