Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

FFT normalization with numpy

Tags:

python

numpy

fft

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!

like image 325
GuntherOnFire Avatar asked Oct 18 '22 02:10

GuntherOnFire


1 Answers

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:

  1. to calculate FFT fft_fwhl = np.fft.fft(fwhl_y)
  2. 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)
  3. 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
  4. 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.

like image 144
GuntherOnFire Avatar answered Oct 21 '22 00:10

GuntherOnFire