how fft points should be rescaled to get the same results as the analytical solution?

Question

I want to work with fast fourier transform using the numpy fft package, and then I am trying to compare the results between the analytical solution and the fast fourier transform, and although I can see with the plots I have done that the curves are similar, it is pretty obvious that the scales are different.

I have tried several different versions of the frequency (angular frequency, frequency and wave number), but all my tries did not work, and in the numpy documentation, it is not clear how the fast fourier transform is exactly defined. For example, I want to work with the fourier transform of the exponential function in time into the angular frequency domain, f(t)=Exp(-a|t|), F(w)=a/pi*(a²+w²) (there are multiple versions of this analytical solution depending on which frequency space we are considering)


def e(t):
    return np.exp(-0.5*abs(t))
def F(w):
    return 0.5/(np.pi)*(1/(((0.5)**2)+((w)**2)))

t=np.linspace(0,100,1000)

w=np.fft.fftfreq(len(t))
plt.plot(w,F(w),'o',label='F(w)')
plt.legend()
plt.show()

fourier=np.fft.fft(e(t))
plt.plot(w,fourier,'o')
plt.show()

I have tried multiple different variants of the above code specially for the frequency, but I am still not getting to the point where the fft and the analytical solution are similar. could anyone please help me out?

francis francis · Accepted Answer · 2019-04-04T19:32:47

The Fourier transform can be applied to integrable functions such as np.exp(-0.5*abs(t)). But the Discrete Fourier Transform computes the Fourier transform of periodic signals. See https://dsp.stackexchange.com/questions/26884/about-fourier-transform-of-periodic-signal and What FFTW Really Computes.

Hence, the DFT of a frame of length T corresponds to the Fourier transform of the periodized frame. Since the frame starts at 0, the Fourier transform of the periodic right sided exponential decay is computed: As you can see, half of the function np.exp(-0.5*abs(t)) is not displayed. Let's correct it and add the periodized increasing part of the two sided exponential decay. I used frequency as a parameter:

import matplotlib.pyplot as plt
import numpy as np

def e(t):
    return np.exp(-0.5*abs(t))
def F(w):
    return 0.5/(np.pi)*(1/(((0.5)**2)+((w)**2)))

def Fc(xi):
    #ok , that's sourced from https://en.wikipedia.org/wiki/Fourier_transform ... Square-integrable functions, one-dimensional, line 207
    return 2*0.5/(((0.5)**2)+(4*(np.pi*xi)**2))


framelength=100.
nbsample=1000
def ep(t):
    #the periodized negative part is added at the end of the frame.
    return np.maximum(np.exp(-0.5*abs(t)),np.exp(-0.5*abs(t-framelength)))

t=np.linspace(0,framelength,nbsample, endpoint=False)

#plotting the periodized signal, to see what's happening
ein=ep(t)
tp=np.linspace(0,10*framelength,10*nbsample, endpoint=False)
periodized=np.zeros(10*nbsample)
for i in range(10):
    for j in range(nbsample):
       periodized[i*nbsample+j]=ein[j]

plt.plot(tp,periodized,'k-',label='periodized frame')
plt.legend()
plt.show()

fourier=np.fft.fft(ep(t))/np.size(ep(t))*framelength

#comparing the mean is useful to check the overall scaling
print np.mean(ep(t))*framelength
print fourier[0]
print Fc(0)

#the frenquencies of the DFT of a frame of length T are 1/T, 2/T ... and negative for the second part of the array.
xi=np.fft.fftfreq(len(t), framelength/len(t))

# comparison between analytical Fourier transform and dft.
plt.plot(xi,Fc(xi),'o',label='F(xi)')
plt.plot(xi,np.real(fourier),'k-', lw=3, color='red', label='DTF')
plt.legend()
plt.show()

And here goes the result:

For an experimental non-periodic signal, an artificial discontinuity appears as the frame is periodized. It induces spectral leakage and windows are applied to attenuate the discontinuity and its effects. One of the potential windows, named a Poisson window, is a two sided exponential decay!

how fft points should be rescaled to get the same results as the analytical solution?

1 Answers