function x = echt(xr, filt_lf, filt_hf, Fs, n)
% ECHT: Endpoint Corrected Hilbert Transform
% X = echt(Xr) computes the so-called discrete-time analytic signal.
% It is a modification of Matlab's hilbert.m function that significantly
% ameliorates the distortions, known as Gibbs phenomenon, at the end of
% the signal.
%
% NOTES:
% 1. One common implementation of the Hilbert Transform uses a DFT
% (aka FFT) as part of its computation. Inherent to the DFT is the
% assumption that a finite sample of a signal is replicated
% infinitely in time, effectively abutting the end of a sample with
% its replicated start. If the start and end of the sample are not
% continuous with each other, distortions are introduced by the DFT.
% Echt effectively smooths out this 'jump discontinuity' by
% selectively deforming the start of the sample. It is hence most
% suited for real-time applications in which the point/s of interest
% is/are the most recent one/s (i.e. last) in the sample window.
% 2. We found that a filter bandwidth (BW= filt_hf- filt_lf) of up to
% half the signal's central frequency works well.
%
% INPUT PARAMETERS:
% xr: time domain signal
% filt_lf: low-cutoff frequency of a bandpass causal filter
% filt_hf: high-cutoff frequency of a bandpass causal filter
% Fs: sampling rate of time domain signal
% n: length of analytic signal (optional)
%
% USE EXAMPLE:
% f0 = 2;
% filt_BW = f0/2;
% N = 1000;
% Fs = N/(2*pi);
% t = -2*pi:1/Fs:0;
% Xr = cos(2*pi*f0*t-pi/4);
% Filt_LF = f0 - filt_BW/2;
% Filt_HF = f0 + filt_BW/2;
% x = echt(Xr, Filt_LF, Filt_HF, Fs, length(Xr))
%
% See also: hilbert, FFT, IFFT, butter, freqz
%
% References:
% Nir Grossman*, David Wang*, Edward S. Boyden "Endpoint Corrected
% Hilbert Transform", IEEE Transaction on Signal Processing (in-review)
% * These authors contributed equally to this work
% Check input
if nargin<5, n=[]; end
if ~isreal(xr)
warning(message('signal:hilbert:Ignore'))
xr = real(xr);
end
% Work along the first nonsingleton dimension
[xr,nshifts] = shiftdim(xr);
if isempty(n)
n = size(xr,1);
end
% Compute FFT
x = fft(xr,n,1); % n-point FFT over columns.
% Set negative components to zero and multiply by 2 postive (apart from DC
% and Nyquist frequency)
h = zeros(n,~isempty(x)); % nx1 for nonempty. 0x0 for empty.
if n > 0 && 2*fix(n/2) == n
% even and nonempty
h([1 n/2+1]) = 1;
h(2:n/2) = 2;
elseif n>0
% odd and nonempty
h(1) = 1;
h(2:(n+1)/2) = 2;
end
x = x.*h(:,ones(1,size(x,2)));
% Compute filter's frequency response
filt_order = 2;
[b,a] = butter(filt_order, [filt_lf filt_hf]/(Fs/2), 'bandpass');
T = 1/Fs*n; % Sampling time period
filt_freq = (ceil(-n/2:n/2-1)'/T); % Center bin frequencies of the FFT.
filt_coeff = freqz(b,a,filt_freq,Fs); % Freq domain filter coefficents.
% Multiply FFT by filter's response function
x = fftshift(x); % rearranges array so DC at the center
x = x.*filt_coeff;
x = ifftshift(x); % rearranges back array so DC at the left
% IFFT
x = ifft(x);
% Convert back to the original shape.
x = shiftdim(x,-nshifts);