HHT1-3要点.docVIP

1、Hilbert边际谱 我觉得既然已经做出EMD了，也就是得到了IMF。这个时候就是做hilbert幅值谱，然后对它积分就可以了。程序不是很难搞到吧！我是用hspec画谱图的，自己又在后面添加了求边际谱的代码，但感觉有问题for k=1:size(E)? ? bjp(k)=sum(E(k,:))*1/fs;??%fs为采样频率;endfigureplot(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: function h1= nspab(data,nyy,minw,maxw,dt) % The function NSPAB generates a smoothed HHT spectrum of data(n,k) % in time-frequency 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.12,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 - 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 specifies 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 Transform data=hilbert(data); a=abs(data); omg=abs(diff(unwrap(a