function [omega,y, xp, yp]=dHif(x);
%
% [omega, y, xp, yp] = dHif(x)
%
% compute the instantaneous frequency of x using the
% discrete Hilbert transform approach documented in 
% D.G. Long "Comments on Hilbert Transform-Based Signal Analysis", MERS
% Report # 04-01, 2004  Available from http://www.mers.byu.edu/
% returns the instanteous frequency estimate omega[n] as well
% as the derivative xp of x, the Hilbert transform y of x and the 
% derivative yp of y.
%
% written 14.2.2004 by DGL at BYU

N=max(size(x));
  
% use short cut computation for faster computation in matlab

oe=x*0;  % odd/even indicator function
der=oe;  % derivative computation function
if N/2==floor(N/2)  % N even
  oe(1:N/2)=1;
  oe(N/2+1:N)=-1;
  der(1:N/2)=0:N/2-1;
  der(N/2+1:N)=(N/2:N-1)-N;
else                % N odd, note center term is left zero
  oe(1:(N-1)/2+1)=1;
  oe((N+1)/2+1:N)=-1;
  der(1:(N-1)/2+1)=0:(N-1)/2;
  der((N+1)/2+1:N)=((N+1)/2:N-1)-N;
end

X=fft(x);
Xr=real(X);Xi=imag(X);

y=real(ifft((Xi-j*Xr).*oe));
xp=real(ifft(j*X.*der))*2*pi/N; %  xp=real(ifft(-(Xi-j*Xr).*der))*2*pi/N;
yp=real(ifft((Xr+j*Xi).*der.*oe))*2*pi/N;

if 1==0    % use explicit formulae to check results
  y1=x*0;
  yp1=y1;
  xp1=y1;
  for n=0:N-1
    sy=0;
    sxp=0;
    syp=0;
    if N/2==floor(N/2)  % N even
      for k=0:N/2-1
	sy=sy+Xr(k+1)*sin(2*pi*n*k/N)+Xi(k+1)*cos(2*pi*n*k/N);
	sxp=sxp-k*Xr(k+1)*sin(2*pi*n*k/N)-k*Xi(k+1)*cos(2*pi*n*k/N);
	syp=syp-k*Xi(k+1)*sin(2*pi*n*k/N)+k*Xr(k+1)*cos(2*pi*n*k/N);
      end  
      for k=N/2:N-1
	sy=sy-Xr(k+1)*sin(2*pi*n*k/N)-Xi(k+1)*cos(2*pi*n*k/N);
	sxp=sxp-(k-N)*Xr(k+1)*sin(2*pi*n*k/N)-(k-N)*Xi(k+1)*cos(2*pi*n*k/N);
	syp=syp+(k-N)*Xi(k+1)*sin(2*pi*n*k/N)-(k-N)*Xr(k+1)*cos(2*pi*n*k/N);
      end
    else   % N odd (center term must be treated as zero)
      for k=0:(N-1)/2
	sy=sy+Xr(k+1)*sin(2*pi*n*k/N)+Xi(k+1)*cos(2*pi*n*k/N);
	sxp=sxp-k*Xr(k+1)*sin(2*pi*n*k/N)-k*Xi(k+1)*cos(2*pi*n*k/N);
	syp=syp-k*Xi(k+1)*sin(2*pi*n*k/N)+k*Xr(k+1)*cos(2*pi*n*k/N);
      end  
      for k=(N+1)/2:N-1
	sy=sy-Xr(k+1)*sin(2*pi*n*k/N)-Xi(k+1)*cos(2*pi*n*k/N);
	sxp=sxp-(k-N)*Xr(k+1)*sin(2*pi*n*k/N)-(k-N)*Xi(k+1)*cos(2*pi*n*k/N);
	syp=syp+(k-N)*Xi(k+1)*sin(2*pi*n*k/N)-(k-N)*Xr(k+1)*cos(2*pi*n*k/N);
      end
    end
    y1(n+1)=sy/N;  % discrete Hilbert transform
    xp1(n+1)=2*pi*sxp/N^2;  % derivative
    yp1(n+1)=2*pi*syp/N^2;  % derivative
  end
  % check explicit vs shortcut computation
  xp_err_max=max(max(abs(xp-xp1)))
  y_err_max=max(max(abs(y-y1)))
  yp_err_max=max(max(abs(yp-yp1)))
end


% compute instantaneous frequency

mag=(x.*x+y.*y);
omega=0*mag;
ind=find(mag ~= 0);
omega(ind)=(yp(ind).*x(ind)-y(ind).*xp(ind))./mag(ind);