EMD分解+三维时频图绘制+伪Wigner-Ville结果对比

1.简单例子

1.1 EMD简单程序

fs = 500;
ts = 1/fs;
t=0:ts:2;
z = sin(2*pi*10*t) + sin(2*pi*5*t)++ sin(2*pi*20*t)++ sin(2*pi*40*t);
imf=emd(z);
emd_visu(z,t,imf)  % EMD专用画图函数

结果  imf = emd(z)  输出六行数据 ，对应 5 个 imf 分量，和 1 个 残差。

f2c 和 c2f 分别从 2 个相反的方向重构原始信号

（ f2c1 是 IMF1 ， f2c2 是 IMF1 和 IMF2 合成的，f2c3 是 IMF1，IMF2，IMF3 合成的，...

    c2f1 是 IMF5 ， c2f2 是 IMF5 和 IMF4 合成的，...）

1.2 三维时频图绘制

%  三维图
[p,q]=ndgrid(t,1:size(imf,1));
plot3(p,q,imf)
grid on
xlabel('Time Values')
ylabel('Mode Number')
zlabel('Mode Amplitude')


%  view(90,0);  % 调整看图角度

 

 

2. EMD相关知识

全称经验模态分解（ Empirical Mode Decomposition,EMD ）,将复杂信号分解为几个具有实际物理意义的固有模态函数（ IMF ），每一个固有模态函数  IMF 均可反映信号本身的内在结构。经验模态分解完全脱离了傅里叶分析的框架，对非平稳态信号处理具有很好的效果。

3.伪Wigner-Ville时频分布（参考help文件）

tfrpwv  伪Wigner-Ville时频分布。

[TFR,T,F] = tfrpwv(X,T,N,H,TRACE)    计算伪Wigner-Ville

离散时间信号X的分布，或交叉两个信号之间的伪Wigner-Ville表示。

X: auto-PWV信号，交叉pwv信号[X1,X2]

T:时间即时( s )          ( 默认值:1:长度(X))。

N:频率箱的个数，默认为长度(X)。

H:频率平滑窗口，在时域内，

H(0)被强制为1(默认值:Hamming(N/4))。

TRACE:如果非零，则显示算法的进程   ( 默认值:0)。

TFR:时频表示。调用时没有输出参数，tfrpwv运行TFRQVIEW。

F:归一化频率向量。

4.伪Wigner-Ville时频结果对比

%计算平滑伪Wigner-Ville分布
sig=hilbert(imf');
sigg=hilbert(z');
% [tfr,t,f]=tfrwv (sig);  %% f 归一化频率向量

subplot(2,1,1);
for i=1:6
  N=length(t);
  [tfr,t,f]=tfrpwv(sig(:,i)); 
  contour(t/fs,f(1:N)*(fs),abs(tfr));
  hold on
end
xlabel('时间 t/s');
ylabel('频率 f/Hz');
title('emd+伪Wigner-Ville')
axis([0 2 0 80])

subplot(2,1,2);
[tfr,t,f]=tfrpwv(sigg);  %% f 归一化频率向量
contour(t/fs,f(1:N)*fs,abs(tfr));  %% 矩阵的等高线图
xlabel('时间 t/s');
ylabel('频率 f/Hz');
title('直接伪Wigner-Ville')
axis([0 2 0 80])