# Approximating Lorentzian Fourier transform with FFT

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')

$$$$


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])
`

• Thanks! I'm interested in recreating both the real and imaginary parts of the lorentzian frequency function separately (i.e real(FTfun), imag(FTfun)), but abs(fft()) takes into account both parts of the fft to create the single lineshape. Is there any way to get both original parts from the fft? – TS1 Aug 9 '19 at 3:24