hspec.m

本人收集的几个版本emd程序

% 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 ) 
%    注意：图形表示是Hilbert能量谱，公式为20 * log ( S * S )
%
% Author: Ivan Magrin-Chagnolleau  <[e-mail]>
% 

function [S,freq] = hspec(imf,N);
%function [S,E,freq] = hspec(imf,N);

%N=500
L = size(imf,1); % Number of components in the decomposition
                 % L 代表分解的基本模式分量数目。
%-------------------------------------------------------------------------
% loop for on each component  对每一个基本模式分量进行循环

S = []; % Matrix which will contain the time-frequency-amplitude representation
        % S 为装载时频幅值表示的矩阵
      
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   解析信号的模在 m 变量中
p = angle(z); % phase of z  解析信号的相位在 p 变量中

for i = 1:L-1   %  attention, here it is justfied by Ren
   
   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(1);
%t=(0:length(x)); % time sample
%fmax=max(max(ceilfreq));
%fmin=min(min(ceilfreq));
%f=t(1:fmax)/length(x)*0.5;% normalized frequency
%f=fliplr(ceilfreq')*10000/N;
%f=(fmin:fmax)*100/N;
  %f=f(0:fmax)*20000/N;     % it is used in yan hua data processing.
         %f=(0:fmax)/N/0.0167;    %used in cracked motor 
%f=f(1:2);
%f=t/length(x)*0.5;  % it cost me too much energe on 0, why?
%     imagesc(t,f,E); % !!! I am not sure it is the best way to visualize it !!!
%contour(t,f,E);
%colorbar;
%set(gca,'YDir','normal');
%xlabel('Time Sample');
%ylabel('Normalized Frequency');
%figure(2);
%freq=1:N/2;
%time=1:length(x)-1;
%mesh(time,freq,E);
      %Marginal=sum(S'.^2);
      %lf=length(f);
   %f2=f(0:fmax)*20000/N;     % it is used in yan hua data processing.
%f2=f(1:lf-1);
%f2=f2+0.006;
%plot(f2,Marginal);
%grid on;
return