function hilbtr_tst(example,l,xl,eps);
%This worksheet studies the Hilbert transform using a discrete
%approximation via the FFT. It provides a collection of functions and their
%Hilbert transforms.
%
clear XL 
%We can study the effects of the length of the sampling interval [-l,l]. If
%the interval is too short, then the approximation is not too good. With
%eps=0, the plot shows the discrete approximation to Hf and the exact
%values obtained by sampling the formula.
%
%We can also study the effect of regularizing the Cauchy principal value
%integral by taking eps to be a small positive number. In this case the
%plot compares the result with the given eps to what is obtained with
%eps=0.
%
%We use 2^k samples
k=16;
if nargin == 2
    XL=l;
    eps=0;
elseif nargin == 3
    eps =0;
    XL=min(l,xl);
elseif nargin == 4
    XL=min(l,xl);
end
%size of the interval over which we sample:
%l=100;
%
%The x-axis limits for the graphs
%xl=10;

%regularization factor. If eps=0, then the graphs of Hf compare to the exact
%result. If eps>0, then the graphs of Hf compare to the result with eps=0.
%eps=0.000;
%Here we choose a function to smooth by choosing example =
%
%example = 5
%
if example == 1
f = @(x) sqrt(1-x.*x).*(abs(x)<=1);
Hfex= @(x) (x+sqrt(x.^2-1)).*(x<-1)+x.*(abs(x)<=1)+(x-sqrt(x.^2-1)).*(x>1);
%
elseif example == 2
f = @(x) real(1./(x+i).^2);
Hfex= @(x) imag(1./(x+i).^2);
elseif example == 3
f = @(x) real(exp(i*x)./(x+i).^2);
Hfex= @(x) imag(exp(i*x)./(x+i).^2);
elseif example == 4
f = @(x) real(1./(1-i*x))
Hfex= @(x) imag(1./(1-i*x));
elseif example == 5
f = @(x) sinc(x);
Hfex= @(x) (1-cos(pi*x))./(pi*x);
end
%Grid for sampling
%
x=linspace(-l,l,2^k);

%
fs=f(x);
figure
subplot(1,2,1)
plot(x,fs)
title(['f'])
Hf=hilbtr(fs,eps);
xlim([-XL,XL])
%
subplot(1,2,2)
if eps == 0
    Hfexs=Hfex(x);
    plot(x,Hf,x,Hfexs)
    xlim([-XL,XL])
else
    Hf0=hilbtr(fs,0);
    plot(x,Hf,x,Hf0)
    xlim([-XL,XL])
end
title('Hf')
end
function Hf=hilbtr(f,eps)
%We use the FFT to compute the Hilbert transform of f which is sampled on 
%a grid 
%
ls=length(f);
ls2=floor((ls+1)/2);
%Compute the FFTs
FFs=fft(f);
%The Fourier multiplier for the Hilbert transform
FHs=[ones(1,ls2) -ones(1,ls-ls2)];
%To study effect of regularization in the convolution integral  we include
%the option to multiply in the frequency domain by exp(-eps*|k|).
reg=exp(-eps*(0:(ls2-1)));
reg=[reg fliplr(exp(-eps*(1:(ls-ls2))))];
FHs=FHs.*reg;
%Multiply and compute the inverse FFT
FHFs=FFs.*FHs;
%This is the discrete approximation to the Hilbert transform:
Hf=real(ifft(FHFs)/i);
end
