Hilbert边际谱.docx

Hilbert边际谱.docx

《Hilbert边际谱.docx》由会员分享，可在线阅读，更多相关《Hilbert边际谱.docx（15页珍藏版）》请在冰豆网上搜索。

Hilbert边际谱

1、Hilbert边际谱

我觉得既然已经做出EMD了，也就是得到了IMF。

这个时候就是做hilbert幅值谱，然后对它积分就可以了。

程序不是很难搞到吧！

我是用hspec画谱图的，自己又在后面添加了求边际谱的代码

fork=1:

size（E）

  bjp（k）=sum（E（k,:

））*1/fs;  %fs为采样频率;

end

figure

plot（bjp）;

xlabel（'频率/Hz'）;

ylabel（'幅值'）;

比如我用两个正弦信号作仿真

fs=1000;

t=1/fs:

1/fs:

1;

y1=2*sin（40*pi*t）;

y2=5*sin（80*pi*t）;

y=[y1,y2];%信号

画出来的图很粗糙，更不用说对实际信号分析了，所以大家看看如何来修正？

黄文章中边际谱对实际信号分析是很好的一条曲线

我用hhspectrum算了一下谱图，同时求了一下边际谱，边际谱程序基本想法同form。

结果也不太好，20HZ处还行，40HZ就有些问题了，见附图

你自己再用这个试试我没有用rilling的hhspectrum

nspab:

functionh1=nspab（data,nyy,minw,maxw,dt）

%ThefunctionNSPABgeneratesasmoothedHHTspectrumofdata（n,k）

%intime-frequencyspace,where

%nspecifiesthelengthoftimeseries,and

%kisthenumberofIMFcomponents.

%Thefrequency-axisrangeisprefixed.

%Negativefrequencysignisreversed.

%

%MATLABLibraryfunctionHILBERTisusedtocalculatetheHilberttransform.

%

%Example,[h,xs,w]=nspab（lod78_p',200,0,0.12,1,3224）.

%

%FunctionsCONTOURorIMGcanbeusedtoviewthespectrum,

%forexamplecontour（xs,w,h）orimg（xs,w,h）.

%

%Callingsequence-

%[h,xs,w]=nspab（data,nyy,minw,maxw,t0,t1）

%

%Input-

%data-2-Dmatrixdata（n,k）ofIMFcomponents

%nyy-thefrequencyresolution

%minw-theminimumfrequency

%maxw-themaximumfrequency

%t0-thestarttime

%t1-theendtime

%Output-

%h-2-DmatrixoftheHHTspectrum,where

%the1stdimensionspecifiesthenumberoffrequencies,

%the2nddimensionspecifiesthenumberoftimevalues

%xs-vectorthatspecifiesthetime-axisvalues

%w-vectorthatspecifiesthefrequency-axisvalues

%Z.Shen（JHU）July2,1995Initial

%-----Getdimensions（numberoftimepointsandcomponents）

[npt,knb]=size（data）;

%-----Gettimeinterval

%-----ApplyHilbertTransform

data=hilbert（data）;

a=abs（data）;

omg=abs（diff（unwrap（angle（data））））/（2*pi*dt）;

%-----Smoothamplitudeandfrequency

filtr=fir1（8,.1）;

fori=1:

knb

a（:

i）=filtfilt（filtr,1,a（:

i））;

omg（:

i）=filtfilt（filtr,1,omg（:

i））;

end

%-----Limitfrequencyandamplitude

fori=1:

knb

fori1=1:

npt-1

ifomg（i1,i）>=maxw,

omg（i1,i）=maxw;

a（i1,i）=0;

elseifomg（i1,i）<=minw,

omg（i1,i）=minw;

a（i1,i）=0;

else

end

end

end

clearfiltrdata

%va=var（omg（200:

1200））

%-----Getlocalfrequency

dw=maxw-minw;

wmx=maxw;

wmn=minw;

%-----Constructtheplotingmatrix

clearp;

h1=zeros（npt-1,nyy+1）;

p=round（nyy*（omg-wmn）/dw）+1;

forj1=1:

npt-1

fori1=1:

knb

ii1=p（j1,i1）;

h1（j1,ii1）=h1（j1,ii1）+a（j1,i1）;

end

end

%-----Do3-pointto1-pointaveraging

[nx,ny]=size（h1）;

%n1=fix（nx/3）;

%h=zeros（n1,ny）;

%fori1=1:

n1

%h（i1,:

）=（h1（3*i1,:

）+h1（3*i1-1,:

）+h1（3*i1-2,:

））;

%end

%clearh1;

%-----Do3-pointssmoothinginx-direction

fltr=1./3*ones（3,1）;

forj1=1:

ny

h1（:

j1）=filtfilt（fltr,1,h1（:

j1））;

end

clearfltr;

%-----Definetheresults

%w=linspace（wmn,wmx,ny-1）';

%xs=linspace（t0,t1,nx）';

h1=flipud（rot90（h1））;

h1=h1（1:

ny-1,:

）;

form求边际谱时所用程序是没有问题的，用的是矩形积分公式。

他所得结果不正确的原因是：

输入的应是调用了toimage后的结果，而不是调用了hhspectrum后的结果。

下面给一段程序，大家可以去试下。

边际谱的分析结果是完全正确的。

clear;

fs=1000;  %fs为采样频率;

N=1000;    %采样点数

t=1/fs:

1/fs:

1;

y1=2*sin（60*pi*t）;

y2=5*sin（90*pi*t）;

y=[y1;y2;zeros（size（y1））];%IMF集

%%%%%%%%%%%%%求边际谱

[A,fa,tt]=hhspectrum（y）;

[E,tt1]=toimage（A,fa,tt,length（tt））;

E=flipud（E）;

fork=1:

size（E,1）

  bjp（k）=sum（E（k,:

））*1/fs;  

end

f=（0:

N-3）/N*（fs/2）;

plot（f,bjp）;

xlabel（'频率/Hz'）;

ylabel（'幅值'）;

%（完整答案）

问：

看了你的程序，我有几点不明白，首先y=[y1;y2;zeros（size（y1））];%IMF集这句代表的含义是什么？

你好像没有作EMD，哪里会有IMF，还有就是E=flipud（E）;，  这句的作用是什么？

答：

一个正弦函数本身就是一个IMF，所以y=[y1;y2;zeros（size（y1））]就是一个IMF集（当然假定了残余函数为0）。

ps：

调用了toimage后得到的结果才是真正的Hilbert谱！

flipud是一个使矩阵上下翻转的函数。

在Grilling提供的程序toimage中，频率是从上往下递增，而通常在时频图中频率应是从下往上递增，所以使用flipud将矩阵翻转后，更便于我们阅读时频图。

对于边际谱来说，如果不对E翻转，边际谱图中的频率将是从从右往左递增的。

2、Hilbert边际谱和FT变换后的幅频谱

这得出的Hilbert边际谱和FT变换后的幅频谱为什么会有这么大的区别呢，到底哪个幅值才是真正的实际幅值呢？

有参考价值吗？

程序如下：

loadshuju

fs=5120;

N=4096;

a1=a（1:

N,1）;

a2=abs（fft（a1））*2/N;

f=fs*（0:

N/2-1）/N;

n=length（f）;

subplot（211）

plot（f,a2（1:

n））

xlabel（'频率/Hz'）;

ylabel（'幅值'）;

title（'FT后的幅频图'）

imf=emd（a1）;

[A,fa,tt]=hhspectrum（imf）;

[E,tt1]=toimage（A,fa,tt,length（tt））;

fork=1:

size（E,1）

  bjp（k）=sum（E（k,:

））*1/fs;

end

f=（0:

N-3）/N*（fs/2）;

subplot（212）

plot（f,bjp）;

xlabel（'频率/Hz'）;

ylabel（'幅值'）;

title（'Hilbert边际谱'）

答：

EMD分解的IMF能量和原信号其能量是不相等的。

所以边际谱能量不能和FFT的能量相比。

边际谱能量只能说明某个信号存在，能量相对于其他的大小。

边际谱是对IMF取包络线。

因此它得到的谱能量要大于被取包络的信号能量。

FFT谱和原信号的能量是相等的。

所以从能量的大小讲，应该是边际谱能量大于FFT的能量。

如果是一个谐波取边际谱和FFT的话，应该是频率对应的能量边际谱大于FFT。

3、EEMD的一些问题

刚刚接触EEMD，从台湾中央大学上下载了程序，然后编写程序如下：

clear;clc

t=1:

1000;

t1=t/100*2*pi;

a1=sin（t1）;

t2=t/10*2*pi;

b1=linspace（0,0,1000）;

fori=250:

350

  b1（i）=0.2*sin（t2（i））;

end

fori=750:

850

  b1（i）=0.2*sin（t2（i））;

end

x=a1+b1;  %x是原信号

subplot（311）;plot（t,a1）;

subplot（312）;plot（t,b1）;

subplot（313）;plot（t,x）;

plot（x）

imf=emd（x）;

emd_visu（x,t,imf）

%eemd

imf_eemd=eemd（x,0.1,100）;

figure

subplot（511）;plot（imf_eemd（:

1））

subplot（512）;plot（imf_eemd（:

2））

subplot（513）;plot（imf_eemd（:

3））

subplot（514）;plot（imf_eemd（:

4））

subplot（515）;plot（imf_eemd（:

5））

figure

subplot（511）;plot（imf_eemd（:

6））

subplot（512）;plot（imf_eemd（:

7））

subplot（513）;plot（imf_eemd（:

8））

subplot（514）;plot（imf_eemd（:

9））

subplot（515）;plot（imf_eemd（:

10））

原信号：

2010-11-1219:

27上传

下载附件（80.08KB）

EMD:

2010-11-1219:

29上传

下载附件（76.54KB）

EEMD:

2010-11-1219:

31上传

下载附件（108.54KB）

2010-11-1219:

31上传

下载附件（76.68KB）

感觉EEMD的第二项为什么会是频率这么高？

和EMD差距这么大？

是不是我程序有问题？

哪位高手请赐教，不胜感激！

答：

我认为第二项应该是随机噪声造成的

虽然集总平均理论上为0

但是毕竟是有限次平均逼近

存在的残余也很正常

我仔细分析了下

图中第一项是信号本身

第二项是残余噪声（随机噪声一般是高频）

第三项是仿真信号中的b1

……

不过这个分解结果于我来看还是不理想

待分解信号中断正弦信号最后没在一个IMF分量中

在没有噪音干扰的情况下，采用EMD得到的分解图可以说是正确的。

我看到一些文献，总是认为EMD分解就是从高频向低频的分解，实际而言，这种理解是对EMD分解得到单分量IMF的一种误解。

如果采用瞬时频率去描述单分量信号，我认为是指信号在某一时刻仅存在一个频率成分。

所以对于原信号而言，[0250]、[350750]与[8501000]区间内的信号只含有10Hz的低频正弦信号，所以在EMD分解过程中，这些信号成分应该与100Hz的高频成分一起出现在第一层IMF中。

第二层IMF中理论上就剩下其余区间内的10Hz正弦信号，但由于EMD分解过程的一些固有缺陷性，才会出现IMF中波形失真等问题。

  对于EEMD，其实际上是利用人为所加的高斯白噪声，来弥补由于频率间断造成的信号间断。

由于高斯白噪声的频率遍布于通频带，这样就相当于原本EMD分解第一层和第二层的那三个区间的正弦信号相对于比其频率更高的噪声频带而言，要远离100Hz的正弦信号，所以从图像上就感觉EEMD解决了所谓的“模态混叠”现象，而实际上与EMD结果区别的部分是被白噪声取代。