% 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 number of frequencies 2^N
df = 1e-2; % frequency step
dt = 1/(Nfft*df); % relation between time and frequency step
fPositive = linspace(df,(Nfft/2)*df,(Nfft/2)); % vector of possitive frequencies (df,fmax)
fData = nan(1,Nfft); % fData is the double sided frequency vector (-fmax,fmax)
fData(Nfft/2+1 - (1:(Nfft/2))) = -fPositive(1:(Nfft/2));
fData((Nfft/2+1)) = 0;
fData((Nfft/2+2):Nfft) = fPositive(1:(Nfft/2-1));
tData = (-Nfft/2*dt:dt:(Nfft/2-1)*dt); % tData is the double sided time vector (-tmax,tmax)

%% These functions form analytical Fourier's pair

H_f = 1./(1 + 1i*(2*pi*fData)); % double sided function of frequency
h_t = exp( - tData).*heaviside(tData); % 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(fData,real(H_f),fData,imag(H_f))
hold on
plot(fData,real(H_f_approx),'x')
hold on
plot(fData,imag(H_f_approx),'x')
grid on

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