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

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))

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))

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函数做出的图较杂乱~