本帖最后由 wdhd 于 2016-3-11 11:14 编辑
不知道为什么会出现错误，请高手指点一下吧
多谢了
- % HSPEC: Hilbert Amplitude Spectrum
- %
- % [S,freq]=hspec(imf,N);
- %
- % S - Time-frequency-amplitude matrix
- % Columns are indexed in time, rows in frequency, values are amplitudes
- %
- % freq- instantaneous frequencies of each component
- %
- % imf - Matrix of intrinsic mode functions (each as a row)
- %
- % N - Number of frequency cells
- %
- % See: Huang et al, Royal Society Proceedings on Math, Physical,
- % and Engineering Sciences, vol. 454, no. 1971, pp. 903-995,
- % 8 March 1998
- %
- % Remark: the graphical representation is the Hilbert Energy Spectrum
- % that is: 20 * log ( S * S )
- %
- % Author: Ivan Magrin-Chagnolleau <ivan@ieee.org>
- %
- function [S,freq] = hspec(imf,N);
- L = size(imf,1); % Number of components in the decomposition
- %-------------------------------------------------------------------------
- % loop for on each component
- S = []; % Matrix which will contain the time-frequency-amplitude representation
- clear x z m p freq
- x = imf'; % now each column is a component
- z = hilbert(x); % analytic signal
- m = abs(z); % module of z
- p = angle(z); % phase of z
- for i = 1:L
- freq(:,i) = instfreq(z(:,i)); % instantaneous frequency
- % if the function instfreq is not available...
- % p(:,i) = unwrap(p(:,i)); % unwrap phase
- % freq(:,i) = abs(diff(p(:,i))); % derivative of the phase and absolute value
- % to have always positive frequencies
- ceilfreq(:,i) = ceil(freq(:,i)*N); % to have integer values - also do a smoothing given the number
- % of frequency cells
- for j = 1:length(x)-2
- S(ceilfreq(j,i),j+1) = m(j+1,i);
- end
- end
- eps = 0.00001; % to avoid zero values before the log
- E = 20 * log ( S.^2 + eps ); % Hilbert energy spectrum
- % plot S
- figure;
- t=1:length(x); % time sample
- f=t/length(x)*0.5; % normalized frequency
- imagesc(t,f,E); % !!! I am not sure it is the best way to visualize it !!!
- colorbar;
- set(gca,'YDir','normal');
- xlabel('Time Sample');
- ylabel('Normalized Frequency');
- return
|
发表于 2006-3-22 11:27
提升卡
变色卡
显身卡
楼主