%    Performs complex signal analysis - estimates instantaneous amplitude, 
%    phase and frequency by taking the Hilbert transform  of a real 
%    band-pass filtered, uniformly sampled time series that is modulating 
%    in frequency and amplitude.
%
%    
%    For real time series s(t) sampled uniformly at dt, the complex signal
%    is S(t)=s(t)+is*(t), where s*(t) = Hilbert transform of s(t).
%    
%    The Hilbert transform is calculated by FFT of time series, zeroing  
%    out Fourier coefficients of negative frequencies, doubling Fourier 
%    coefficients of positive frequencies, followed by inverse FFT.
%    
%    Instantaneous amplitude: iamp(t) = |s(t)+is*(t)|
%    Instantaneous phase:     ifaze(t) = unwrap(arctan[s*(t)/s(t)])
%    Instantaneous frequency: ifreq(t) = diff(ifaze(t))/(2*pi*dt)
%
%    INPUTS:
%           s=[t,x] -- bandpassed, centered, uniformly sampled time series 
%                      with time column and value column
%           
%    OUTPUTS:
%           iamp = instantaneous amplitude (in units of x)
%           ifaze = instantaneous phase (in radians)
%           t = time column of input time series
%           ifreq = instantaneous frequency (in cycles/(time unit))
%           tm1=time column reduced by 1 for ifreq
%
%      2002: Designed by Florian Maurer.
%      2006: Updated by Linda Hinnov.  
%      2017: Rewritten with explicit Hilbert Transform by Linda Hinnov.
%
%      References:
%        
%   Kodama, K.P. and Hinnov, L.A. (2015), Chapter 4 in Rock Magnetic
%        Cyclostratigraphy, Wiley and Sons.
%   Taner, M.T., Koehler, F., & Sheriff, R.E. (1979) Complex trace 
%   analysis. Geophysics,44, 1041?1063. DOI:10.1190/1.1440994.
%
%   %TEST:
%
% Load "hilbertsignaltest.mat", note the variable "signal" in Workspace, 
% and in the Command Window enter:
%
% [iamp,ifaze,t,ifreq,tm1]=hilbertsignal(signal);
%
function [iamp,ifaze,t,ifreq,tm1]=hilbertsignal(s)
npts = length(s);
t=s(:,1);
dt=abs(t(2)-t(1));
a=s(:,2);
b=zeros(length(a),1);
dt=abs(t(2)-t(1));
% Store real signal in complex variable
scom=complex(a,b);
% Forward FFT; honor time series length npts
scomfft=fft(scom,npts);
% Divide into real a and imaginary b
a=real(scomfft);
b=imag(scomfft);
% No. positive frequencies odd or even
npts2=fix(npts/2);
xnpts=npts;
xnpts2=xnpts/2;
diff=0.0;
diff = xnpts2-npts2;
if diff~=0.0 
npts2=npts2+1;
end;
% Double the positive frequency values
for n=1:npts2
a(n)=2*a(n);
b(n)=2*b(n);
end;
% Zero out the negative frequency values
npts2p1=npts2+1;
for n=npts2p1:npts
a(n)=0.0;
b(n)=0.0;
end
% Inverse FFT, honor input time series length npts
scomfft=complex(a,b);
scom=ifft(scomfft,npts);
% Calculate instantaneous amplitude
a=real(scom);
b=imag(scom);
iamp = a.^2 + b.^2;
iamp = sqrt(iamp);
% Calculate instantaneous phase, unwrap and detrend
faze=angle(scom);
ufaze=unwrap(faze);
length(ufaze);
ifaze=detrend(ufaze);
% Calculate instantaneous frequency
ifreq=[ ];
npts1=npts-1;
for n=1:npts1
n1=n+1;
ifreq(n)=(ufaze(n1)-ufaze(n))/(2*pi*dt);
end
length(ifreq);
tm1=t(1:(length(t)-1)); % time column for ifreq; one less due to diff
% Plot the results
figure;
subplot(4,1,1), plot(t,s(:,2)),title('Input modulating signal')
subplot(4,1,2), plot(t,iamp),title('Instantaneous amplitude')
subplot(4,1,3), plot(t,ifaze),title('Instantaneous phase')
subplot(4,1,4), plot(tm1,ifreq),title('Instantaneous frequency')
% EOF



