function [z,A,phase] = hilbert(x)
% Hilbert transform
% Routine for computing Hilbert transform according
% to Bendat and Piersol.
% z = analytic signal = x + i*HT(x)
% A = envelope function
% phase = instantaneous phase
%
% by D. Holmgren, August 16 94 holmgren@phobos.astro.uwo.ca
X = fft(x);
nl = length(x); no2 = nl/2; n1 = nl - 1;
Z = zeros(nl,1);
Z(1:no2) = X(1:no2);
% construct analytic signal, imaginary part of which is
% the desired Hilbert transform...
z = 2 .* ifft(Z);
z = fftshift(z);
% compute amplitude and phase of z...
A = sqrt(real(z).^2 + imag(z).^2);
phase = atan(imag(z) ./ real(z));