1

I have a temporal signal and I calculate its Fourier Transform to get the frequencial signal. According to Parseval's theorem, the two signals have the same energy. I successfully demonstrate it with Python. However, when I calculate the inverse Fourier Transform of the frequencial signal, the energy is no longer conserved. Here is my code:

import numpy as np
import numpy.fft as nf
import matplotlib.pyplot as plt

#create a gaussian as a temporal signal    
x = np.linspace(-10.0,10.0,num=1000)
dx = x[1]-x[0]
sigma = 0.4
gx = (1.0/(2.0*np.pi*sigma**2.0)**0.5)*np.exp(-0.5*(x/sigma)**2.0)

#calculate the spacing of the frequencial signal
f=nf.fftshift(nf.fftfreq(1000,dx))
kk = f*(2.0*np.pi)
dk = kk[1]-kk[0]

#calculate the frequencial signal (FT)
#the convention used here allows to find the same energy
gkk = nf.fftshift(nf.fft(nf.fftshift(gx)))*(dx/(2.0*np.pi)**0.5)

#inverse FT
gx_ = nf.ifft(nf.ifftshift(gkk))*dk/(2 * np.pi)**0.5

#Parseval's theorem
print("Total energy in time domain = "+str(sum(abs(gx)**2.0)*dx))
print("Total energy in freq domain = "+str(sum(abs(gkk)**2.0)*dk))
print("Total energy after iFT = "+str(sum(abs(gx_)**2.0)*dx))

After executing this code, you can see that the two first energies are the same, whereas the third is orders magnitude less than the two first, although I am supposed to find the same energy. What happened here?

gravedigger
  • 359
  • 1
  • 3
  • 13

1 Answers1

5

The numpy FFT procedures actually and in contrast to other software do adjust for the sequence length, so that you get

nf.ifft(nf.fft(gx)) == gx

up to some floating point error. If your dx and dk are computed the usual way, then dk*dx=(2*pi)/N which only works for unadjusted FFT routines.

You can test the behavior of numpy.fft using

In [20]: sum(abs(gx)**2.0)
Out[20]: 35.226587122763036

In [21]: gk = nf.fft(gx)

In [22]: sum(abs(gk)**2.0)
Out[22]: 35226.587122763049

In [23]: sum(abs(nf.ifft(gk))**2.0)
Out[23]: 35.226587122763014

which tells us that the fft is the usual unadjusted transform and ifft divides the result by sequence length N=num. The typical ifft can be emulated by

gxx = (nf.fft(gk.conj())).conj()

then you get that

gx == gxx/1000

up to floating point errors. Or you can reverse the adjustment using

#inverse FT
gx_ = nf.ifft(nf.ifftshift(gkk))*(num*dk)/(2 * np.pi)**0.5
Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • Thank you for your answer and the detailed explanation! Sorry for the incomplete code that I included in my question. Where this adjustment comes from? Why there is no such an adjustment in continuous inverse FT formula? – gravedigger May 06 '15 at 11:53
  • You get the constant by studying the transform of the constant signal with value 1. In the continuous case, this gives you the infinity of the delta peak. – Lutz Lehmann May 06 '15 at 12:00