HHT13讲解Word格式文档下载.docx

1、y=y1,y2; % 信号画出来的图很粗糙，更不用说对实际信号分析了，所以大家看看如何来修正？黄文章中边际谱对实际信号分析是很好的一条曲线我用hhspectrum算了一下谱图，同时求了一下边际谱，边际谱程序基本想法同form。结果也不太好，20HZ处还行，40HZ就有些问题了，见附图你自己再用这个试试 我没有用rilling的hhspectrumnspab:function h1= nspab（data,nyy,minw,maxw,dt）% The function NSPAB generates a smoothed HHT spectrum of data（n,k） % in time-f

2、requency space, where % n specifies the length of time series, and % k is the number of IMF components.% The frequency-axis range is prefixed.% Negative frequency sign is reversed.% MATLAB Library function HILBERT is used to calculate the Hilbert transform.% Example, h,xs,w = nspab（lod78_p,200,0,0.1

3、2,1,3224）.% Functions CONTOUR or IMG can be used to view the spectrum,% for example contour（xs,w,h） or img（xs,w,h）.% Calling sequence-% h,xs,w = nspab（data,nyy,minw,maxw,t0,t1）% Input-% data - 2-D matrix data（n,k） of IMF components% nyy - the frequency resolution% minw - the minimum frequency% maxw

4、- the maximum frequency% t0 - the start time% t1 - the end time% Output-% h - 2-D matrix of the HHT spectrum, where% the 1st dimension specifies the number of frequencies,% the 2nd dimension specifies the number of time values% xs - vector that specifies the time-axis values% w - vector that specifi

5、es the frequency-axis values% Z. Shen （JHU） July 2, 1995 Initial%- Get dimensions （number of time points and components）npt,knb = size（data）;%- Get time interval%- Apply Hilbert Transformdata=hilbert（data）;a=abs（data）;omg=abs（diff（unwrap（angle（data）/（2*pi*dt）;%- Smooth amplitude and frequency filtr=

6、fir1（8,.1）; for i=1:knb a（:,i）=filtfilt（filtr,1,a（:,i）; omg（:,i）=filtfilt（filtr,1,omg（:%- Limit frequency and amplitudeknb for i1=1:npt-1 if omg（i1,i） =maxw, omg（i1,i）=maxw; a（i1,i）=0; elseif omg（i1,i）=minw, omg（i1,i）=minw; else endclear filtr data%va=var（omg（200:1200）%- Get local frequencydw=maxw -

7、 minw;wmx=maxw;wmn=minw;%- Construct the ploting matrixclear p;h1=zeros（npt-1,nyy+1）;p=round（nyy*（omg-wmn）/dw）+1;for j1=1: ii1=p（j1,i1）; h1（j1,ii1）=h1（j1,ii1）+a（j1,i1）;%- Do 3-point to 1-point averagingnx,ny=size（h1）;%n1=fix（nx/3）;%h=zeros（n1,ny）;%for i1=1:n1 %h（i1,:）=（h1（3*i1,:）+h1（3*i1-1,:）+h1（3*i

8、1-2,:）;%end%clear h1;%- Do 3-points smoothing in x-directionfltr=1./3*ones（3,1）;ny h1（:,j1）=filtfilt（fltr,1,h1（:,j1）;clear fltr;%- Define the results%w=linspace（wmn,wmx,ny-1）;%xs=linspace（t0,t1,nx）h1=flipud（rot90（h1）;h1=h1（1:ny-1,:form求边际谱时所用程序是没有问题的，用的是矩形积分公式。他所得结果不正确的原因是：输入的应是调用了toimage后的结果，而不是调用了

9、hhspectrum后的结果。下面给一段程序，大家可以去试下。边际谱的分析结果是完全正确的。clear; %fs为采样频率;N=1000;%采样点数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）;size（E,1）f=（0:N-3）/N*（fs/2）;plot（f,bjp）;%（完整答案）问：看了你的程序，我有几点不明白，首先y=y1; %IMF集这句代表的含义是什么？你

10、好像没有作EMD，哪里会有IMF，还有就是E=flipud（E）; ，这句的作用是什么？答：一个正弦函数本身就是一个IMF，所以y=y1;zeros（size（y1）就是一个IMF集（当然假定了残余函数为0）。ps：调用了toimage后得到的结果才是真正的Hilbert谱！flipud是一个使矩阵上下翻转的函数。在Grilling 提供的程序toimage中，频率是从上往下递增，而通常在时频图中频率应是从下往上递增，所以使用flipud将矩阵翻转后，更便于我们阅读时频图。对于边际谱来说，如果不对E翻转，边际谱图中的频率将是从从右往左递增的。2、Hilbert边际谱和FT变换后的幅频谱这得出的

11、Hilbert边际谱和FT变换后的幅频谱为什么会有这么大的区别呢，到底哪个幅值才是真正的实际幅值呢？有参考价值吗？ 程序如下：load shujufs=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）title（FT后的幅频图）imf=emd（a1）;A,fa,tt=hhspectrum（imf）;subplot（212）Hilbert边际谱EMD分解的IMF能量和原信号其能量是不相等的。所以边际谱能量不能和FFT的能量相比。边际谱能量只能说明某个

12、信号存在，能量相对于其他的大小。边际谱是对IMF取包络线。因此它得到的谱能量要大于被取包络的信号能量。FFT谱和原信号的能量是相等的。所以从能量的大小讲，应该是边际谱能量大于FFT的能量。如果是一个谐波取边际谱和FFT的话，应该是频率对应的能量边际谱大于FFT。3、EEMD的一些问题刚刚接触EEMD，从台湾中央大学上下载了程序，然后编写程序如下：clct=1:1000;t1=t/100*2*pi;a1=sin（t1）;t2=t/10*2*pi;b1=linspace（0,0,1000）;for i=250:350 b1（i）=0.2*sin（t2（i）;for i=750:850x=a1+b1

13、;%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）%eemdimf_eemd=eemd（x,0.1,100）;subplot（511）;plot（imf_eemd（:,1）subplot（512）;,2）subplot（513）;,3）subplot（514）;,4）subplot（515）;,5）,6）,7）,8）,9）,10）原信号：2010-11-12 19:27 上传下载附件 （80.08 KB） EMD:29

14、上传下载附件 （76.54 KB） EEMD:31 上传下载附件 （108.54 KB） 下载附件 （76.68 KB） 感觉EEMD的第二项为什么会是频率这么高？和EMD差距这么大？是不是我程序有问题？哪位高手请赐教，不胜感激！我认为第二项应该是随机噪声造成的虽然集总平均理论上为0但是毕竟是有限次平均逼近存在的残余也很正常我仔细分析了下图中第一项是信号本身第二项是残余噪声（随机噪声一般是高频）第三项是仿真信号中的b1不过这个分解结果于我来看还是不理想待分解信号中断正弦信号最后没在一个IMF分量中在没有噪音干扰的情况下，采用EMD得到的分解图可以说是正确的。我看到一些文献，总是认为EMD分解就

15、是从高频向低频的分解，实际而言，这种理解是对EMD分解得到单分量IMF的一种误解。如果采用瞬时频率去描述单分量信号，我认为是指信号在某一时刻仅存在一个频率成分。所以对于原信号而言，0 250、350 750与850 1000区间内的信号只含有10Hz的低频正弦信号，所以在EMD分解过程中，这些信号成分应该与100Hz的高频成分一起出现在第一层IMF中。第二层IMF中理论上就剩下其余区间内的10Hz正弦信号，但由于EMD分解过程的一些固有缺陷性，才会出现IMF中波形失真等问题。 对于EEMD，其实际上是利用人为所加的高斯白噪声，来弥补由于频率间断造成的信号间断。由于高斯白噪声的频率遍布于通频带，这样就相当于原本EMD分解第一层和第二层的那三个区间的正弦信号相对于比其频率更高的噪声频带而言，要远离100Hz的正弦信号，所以从图像上就感觉EEMD解决了所谓的“模态混叠”现象，而实际上与EMD结果区别的部分是被白噪声取代。