fft - Proving Fourier transform operation in Python

Question

I have an expression in the time domain

f = -1j*H(t) * exp(-(1j*a+b)*t)

which can be Fourier transformed analytically using known properties (H is the Heaviside step function). The result of this FT operation is

F = (w-a-1j*b)/((w-a)**2+b**2)

where w is frequency.

Now I'm using the tips in this article to do numerical Fourier transform on f in Python, and confirm that I do get the same analytical result F:

import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(-10,10,1e4) # time
w = np.linspace(-10,10,1e4) # frequency

b = 0.1
a = 1

H = lambda x: 1*(x>0) # heaviside function

# function in time
f = -1j*H(t)*np.exp(-(1j*a+b)*t)

# function in frequency (analytical work)
F = (w-a-1j*b)/((w-a)**2+b**2)

hann = np.hanning(len(t))  # hanning window

# function in frequency (numerical work)
F2 = 2/len(t)*np.fft.fft(hann*f)

plt.figure()
plt.plot(w,F.real,'b',label='analytical')
plt.plot(w,F2.real,'b--',label='fft')
plt.xlabel(r'$\omega$')
plt.ylabel(r'Re($F(\omega)$)')
plt.legend(loc='best')

plt.figure()
plt.plot(w,F.imag,'g',label='analytical')
plt.plot(w,F2.imag,'g--',label='fft')
plt.xlabel(r'$\omega$')
plt.ylabel(r'Im($F(\omega)$)')
plt.legend(loc='best')

plt.show()

However Python's FFT function seems to give me something completely wrong. This is evident when F and F2 are plotted.

Edit: Here are the plots...

It's not obvious in these figures, but if you zoom in around the w=-10 and 10 areas, there are small oscillations, possibly due to the fft algorithm.

Cris Luengo Cris Luengo · Accepted Answer · 2018-03-07T01:49:05

The FFT algorithm computes the DFT, which has the origin (both spatial and in frequency domain) on the first sample. You need to shift your signal (after applying the Hanning window) so that t=0 is the leftmost sample, and after computing the FFT you have to do the inverse shift.

MATLAB has ifftshift and fftshift, which implement those two shifts. NumPy must have similar functions.

Another issue with your code is that you compute the DFT, and plot it at the locations given by the w that you computed, but is unrelated to the actual frequencies at which the DFT is computed.

Here is your code, translated to MATLAB, and fixed to properly compute F2 and w *. I hope this is useful. One thing to note is that your F does not match F2, I am confident that this is not due to an error in F2, but an error in your computation of F. The shapes are similar, but F is scaled differently and mirrored.

N = 1e3;
t = linspace(-100,100,N); % time
Fs = 1/(t(2)-t(1));
w = Fs * (-floor(N/2):floor((N-1)/2)) / N; % NOTE proper frequencies

b = 0.1;
a = 1;

H = @(x)1*(x>0); % Heaviside function

% function in time
f = -1j*H(t).*exp(-(1j*a+b)*t);

% function in frequency (analytical work)
F = (w-a-1j*b)./((w-a).^2+b.^2);

% hanning window
hann = 0.5*(1-cos(2*pi*linspace(0,1,N)));

% function in frequency (numerical work)
F2 = fftshift(fft(ifftshift(hann.*f))); % NOTE shifting of origin

figure

subplot(2,1,1), hold on
plot(w,real(F),'b-')
plot(w,real(F2),'r-')
xlabel('\omega')
ylabel('Re(F(\omega))')
legend({'analytical','fft'},'Location','best')

subplot(2,1,2), hold on
plot(w,imag(F),'b-')
plot(w,imag(F2),'r-')
xlabel('\omega')
ylabel('Im(F(\omega))')
legend({'analytical','fft'},'Location','best')

^{Footnote:

* I also changed the colors, MATLAB's green is too light.}

fft - Proving Fourier transform operation in Python

1 Answers