MATLAB实现EMD分解及希尔伯特谱分析

MATLAB实现EMD分解及希尔伯特谱分析

希尔伯特—黄变换

传统的傅里叶变化只能得到信号在采样周期内的全局频率信息，处理频率随时间变化的非平稳信号具有很大的局限性，希尔伯特－黄变换(Hilbert-Huang Trans form) 是由N. E. Huang 等人于1998年提出的一种非线性、非平稳信号的分析处理方法，可以求得信号的瞬时频率、瞬时幅值等瞬时特征。

原理

希尔伯特变换：

$\stackrel{\wedge }{g}(t)=\frac{1}{\pi }P\underset{-\infty }{\overset{+\infty }{\int }}\frac{g(\tau )}{t-\tau }d\tau =\frac{P}{\pi t}\ast g(t)$

信号g(t) 的希尔伯特变换是原信号g(t) 与1/$\pi 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);
imf=imf';
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函数做出的图较杂乱~