% [A,f,tt] = HHSPECTRUM2(imf) computes the Hilbert-Huang spectrum
%
% inputs:
% 	- imf : matrix with one IMF per row
%
% outputs:
%   - A   : amplitudes of IMF rows
%   - f   : instantaneous frequencies
%   - tt  : truncated time instants
%
% calls:
%   - hilbert  : computes the analytic signal
%

function [A,f,tt] = hhspectrum2(imf)

t=1:size(imf,2);
lt = length(t);
tt = t(1:lt-2);

N=128;

for i=1:(size(imf,1)-1)
  
  x = imf(i,:)';
  
  % blackman on both ends
  w = blackman(2*N);
  x(1:N) = x(1:N) .* w(1:N);
  x(end-N+1:end) = x(end-N+1:end) .* w(end-N+1:end);
 
  % add zeros
  x = [x; zeros(size(x))];

  % Hilbert
  xa = hilbert(x);

  % remove zeros
  xa(length(xa)/2+1:end) = [];
  
  an(i,:) = xa;
  
  t=2:length(xa)-1;
  
  % diff of angles
  fnormhat= angle(-xa(t+1).*conj(xa(t-1)))+pi;

  % amplitude
  A(i,:)=abs(an(i,2:end-1));

  % average
  a2 = A(i,:) .* A(i,:);
  num = conv(w, fnormhat' .* a2);
  den = conv(w, a2);
  
  % frequencies
  f(i,:) = 0.5*wkeep((num ./ den), lt-2)/(2*pi);

end
