function [A,f,tt] = hhspectrum(imf,t,l,aff)% [A,f,tt] = HHSPECTRUM(imf,t,l,aff) computes the Hilbert-Huang spectrum%% inputs:% 	- imf : matrix with one IMF per row%   - t   : time instants%   - l   : estimation parameter for instfreq%   - aff : if 1, displays the computation evolution%% outputs:%   - A   : amplitudes of IMF rows%   - f   : instantaneous frequencies%   - tt  : truncated time instants%% calls:%   - hilbert  : computes the analytic signal%   - instfreq : computes the instantaneous frequencyif nargin < 2  t=1:size(imf,2);endif nargin < 3  l=1;endif nargin < 4  aff = 0;endlt=length(t);tt=t((l+1):(lt-l))';for i=1:(size(imf,1)-1)  an(i,:)=hilbert(imf(i,:)')';% %---------------------------------% x=an(i,:)';% theta=angle(x);% % omega=diff(theta); % instantanous frequency% L_omega=size(omega);% % for k=1:L_omega%     if omega(k)<=0%         omega(k)=omega(k)+2*pi;   %if instantanous frequency less than zero, add 2*pi to omega%     else%     end% end% fff=0.5*omega/(2*pi);% for j=2:lt%     fnormhat(i,j-1)=(0.5*(angle(-x(j).*conj(x(j-1)))+pi)/(2*pi))';% end% %---------------------------------   f(i,:)=instfreq(an(i,:)')';      A=abs(an(:,l+1:end-l));  if aff    disp(['mode ',int2str(i),'trait'])  endend