Approximating Lorentzian Fourier transform with FFT

Question

I have a Lorentzian frequency distribution

$F(w) = \frac{1+iz}{1+z^2}$

Where

$z = \frac{w-\Omega}{R}$

With $\Omega$ being the peak frequency and R the decay constant. I know that analytically the Fourier transform should be

$F(t) = exp(i\Omega 2\pi t)exp(-Rt)$

When i take the FFT of this expression, it doesn't return the original frequency plot. I understand that there may be a scaling factor (1/n) that should be in there somewhere, but even when i scale for one frequency, if i then change $\Omega$ or R, the amplitude is no longer scaled properly, suggesting that the scaling factor is a function of $\Omega$ and/or R. The FFT also seems to be mirrored along the intensity axis.

I'm fairly new to DSP, but i do understand that the continuous fourier transform is not the discrete fourier transform. I've read this (https://dspillustrations.com/pages/posts/misc/approximating-the-fourier-transform-with-dft.html) but that approach makes the approximation worse.

I'd like to return the original frequency distribution when i take the FFT of my time signal. Am i missing something fundamental or is it a fairly simple scaling error? I've attached my code below.

Cheers.

# R script to compare FFT and Analytical fourier transform
library(SynchWave)

#-------------------------------------------------
# Frequency and time axes
n <- 100
f <- seq(0, 1, length.out = n)
t <- seq(0, n, length.out = n)

# peak paramaters 
O <- 0.3 # Frequency values from 0->1
R <- 0.04 # Decay in arbritrary units

z <- (f-O)/R

# The original lorentzian frequency 
ff <-complex(re = 1, im = z)/(1 + z^2)

# creating the time domain signal
ftideal <- exp(-R*t)*exp(complex(i = (O)*2*pi*t))

unscaled <- (fft(ftideal))
scaled <- unscaled - min(Re(unscaled))

plot(f, Re(ff), type = 'l')
lines(f, Re(scaled), type = "l", col = 'red')

```

Answer 1

First of all, take a look at this very related question and its answer(s). Second, your analytic solution is wrong, so it's no surprise that you don't see the correct result.

The time domain function

$$x(t)=e^{j2\pi f_0t}e^{-Rt}\tag{1}$$

has the following Fourier transform:

$$X(f)=\frac{1}{R}\frac{1-j\frac{2\pi(f-f_0)}{R}}{1+\left(\frac{2\pi(f-f_0)}{R}\right)^2}\tag{2}$$

Also, scaling should be multiplicative and not additive (unless you're in the logarithmic domain), and you should also not only look at the real parts of complex-valued functions.

As explained in this answer, approximating the CTFT by the DFT usually results in two types of errors: the truncation error (due to truncation of the time domain function), and the aliasing error (due to sampling the time domain function). These errors can be made small by choosing a large sampling frequency and a large DFT length.

Below is an example code in Octave/Matlab showing how the CTFT of the given function can be approximated by the DFT:

F0 = 0.3;
R = 0.04;

Fs = 30;            % sampling frequency
Ts = 1/Fs;
Tmax = 100;         % length of time domain signal
N = round(Tmax/Ts); % DFT length

t = (0:N-1) * Tmax / (N-1);
N2 = round(N/2);
f = (0:N2-1)/N*Fs;

fun = exp(1i*2*pi*F0*t) .* exp(-R*t);

% analytic expression for CTFT
z = 2*pi*(f - F0) / R;
FTfun = (1 - 1i*z) ./ (1 + z.^2) / R;

% DFT approximation of CTFT
FTfun2 = Ts * fft(fun,N);
FTfun2 = FTfun2(1:N2);

plot(f,20*log10(abs(FTfun)),f,20*log10(abs(FTfun2)))
    title('magnitudes (dB)'), xlabel('f'), legend('CTFT','DFT'), grid on
    axis([0,Fs/2,-40,30])