FFT normalization with numpy

Question

Just started working with numpy package and started it with the simple task to compute the FFT of the input signal. Here's the code:

import numpy as np
import matplotlib.pyplot as plt

#Some constants
L = 128
p = 2
X = 20
x = np.arange(-X/2,X/2,X/L)
fft_x = np.linspace(0,128,128, True)

fwhl = 1

fwhl_y = (2/fwhl) \
*(np.log([2])/np.pi)**0.5*np.e**(-(4*np.log([2]) \
*x**2)/fwhl**2)

fft_fwhl = np.fft.fft(fwhl_y, norm='ortho')

ampl_fft_fwhl = np.abs(fft_fwhl)

plt.bar(fft_x, ampl_fft_fwhl, width=.7, color='b')

plt.show()

Since I work with an exponential function with some constant divided by pi before it, I expect to get the exponential function in Fourier space, where the constant part of the FFT is always equal to 1 (zero frequency). But the value of that component I get using numpy is larger (it's about 1,13). Here I have an amplitude spectrum which is normalized by 1/(number_of_counts)**0.5 (that's what I read in numpy documentation). I can't understand what's wrong... Can anybody help me?

Thanks!

[EDITED] It seems like the problem is solved, all you need to get the same result of Fourier integral and of FFT is to multiply FFT by the step (in my case it's X/L). And as for normalization as option of numpy.fft.fft(..., norm='ortho'), it's used only to save the scale of the transform, otherwise you'll need to divide the result of the inverse FFT by the number of samples. Thanks everyone for their help!

Be careful with division - if you divide two integers in python they don't return a float. If you divide 15 by 10, you'll get 1, not 1.5. Similarly, 20/128= 0, not .15 — flyingmeatball
@flyingmeatball That's true on python 2.x. On python 3.x they are casted to float explicitely ;) — BPL
The zero frequency corresponds to the mean of the input: fft_fwhl[0] # (0.56568542494923801+0j) which matches np.mean(fwhl_y) * np.sqrt(len(fwhl_y)) # 0.56568542494923812 (it’s sqrt because you’re using norm="ortho"). — Ahmed Fasih
I find it more helpful to think of FFT in terms of the DFT definition, instead of the continuous Fourier transform: docs.scipy.org/doc/numpy/reference/… it’s often a little difficult to get the DFT/FFT to give the same answer as the continuous Fourier transform. — Ahmed Fasih
For everyone looking for a more thorough introduction to FFT scaling (and windowing!), I can recommend this paper by Heinzel et al. which helped me a lot to understand FFT scaling. It's a lot more application oriented than most introductory texts and still has enough theory to go beyond a simple recipe: Spectrum and spectral density estimation by the Discrete Fouriertransform (DFT), including a comprehensive list of window functions and some new flat-top windows, G. Heinzel et al., 2002 — cx05

GuntherOnFire GuntherOnFire · Accepted Answer · 2016-09-06T10:05:01

I've finally solved my problem. All you need to bond FFT with Fourier integral is to multiply the result of the transform (FFT) by the step (X/L in my case, FFTX/L), it works in general. In my case it's a bit more complex since I have an extra rule for the function to be transformed. I have to be sure that the area under the curve is equal to 1, because it's a model of δ function, so since the step is unchangeable, I have to fulfill stepsum(fwhl_y)=1 condition, that is X/L=1/sum(fwhl_y). So to get the correct result I have to make following things:

to calculate FFT fft_fwhl = np.fft.fft(fwhl_y)
to get rid of phase component which comes due to the symmetry of fwhl_y function, that is the function defined in [-T/2,T/2] interval, where T is period and np.fft.fft operation thinks that my function is defined in [0,T] interval. So to get amplitude spectrum only (that's what I need) I simply use np.abs(FFT)
to get the values I expect I should multiply the result I got on previous step by X/L, that is np.abs(FFT)*X/L
I have an extra condition on the area under the curve, so it's X/L*sum(fwhl_y)=1 and I finally come to np.abs(FFT)*X/L = np.abs(FFT)/sum(fwhl_y)

Hope it'll help anyone at least.

FFT normalization with numpy

2 Answers