% Numerical evaluation of Fourier's transform via FFT
% h(t) <---> H(f), 2*pi*f=omega, (FFT is natively made in f, not omega)

clear all
close all
clc

Nfft=2^15; % Nfft is the number of frequencies, 2^N
df=1e-2; % frequency step
dt=1/(Nfft*df); % relation between time and frequency step
f_positive=linspace(df,(Nfft/2)*df,(Nfft/2)); % vector of possitive frequencies (df,fmax)
f_data=nan(1,Nfft); % f_data is a double sided frequency vector (-fmax,fmax)
f_data(Nfft/2+1-(1:(Nfft/2)))=-f_positive(1:(Nfft/2));
f_data((Nfft/2+1))=0;
f_data((Nfft/2+2):Nfft)=f_positive(1:(Nfft/2-1));
t_data=(-Nfft/2*dt:dt:(Nfft/2-1)*dt); % t_data is a double sided time vector (-tmax,tmax)

%% These functions form an analytical Fourier's pair

H_f=1./(1+1i*(2*pi*f_data)); % double sided function of frequency
h_t=exp(-t_data).*heaviside(t_data); % double sided function of time

%% Numerical evaluation of Fourier's transform

h_t_approx =fftshift( ifft( ifftshift( H_f ) ) ) / dt; % inverse Fourier transform
H_f_approx =fftshift( fft( ifftshift( h_t ) ) ) * dt; % Fourier transform

%% Comparison

figure
plot(f_data,real(H_f),f_data,imag(H_f))
hold on
plot(f_data,real(H_f_approx),'-x')
hold on
plot(f_data,imag(H_f_approx),'-x')
grid on

figure
plot(t_data,real(h_t),t_data,imag(h_t))
hold on
plot(t_data,real(h_t_approx),'-x')
hold on
plot(t_data,imag(h_t_approx),'-x')
grid on