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)=π1P−∞∫+∞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:Nfor j=1:N-1if((df_c*i)>df(j)&&(df_c*(i-1))<df(j))C(i,j)=E(j);endend
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:3subplot(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:Nfor j=1:N-1if((df_c*i)>df(j)&&(df_c*(i-1))<df(j))C(i,j)=E(j);endendend
endfigure(4)
imagesc([0 tmax],[0 fs/10],C)
xlabel('t');
ylabel('f');
title('hht')
% imagesc(tt,df,E)
colorbarfigure(5)
hht(imf,fs,'FrequencyLimits',[0 50])

经对比可以发现，自创函数做出的图（需要倒过来一下）与matlab自带的hht函数相比，matlab自带的hht函数做出的图较杂乱~