pythonnumpyfftinverse-transform

inverse Fourier unclear Data


I am currently trying to calculate the inverse Fourier of a known function. Unfortunately the np.fft.ifft() function did not show the result I was looking for. Since I was not sure where the error accured I constructed a known problem and was confronted with a similar error. The Fourier function was 1/(1+x^2) and the inverse Fourier of this function is a*exp(-|k|). With a being a constant. As seen in the Plot the data produced via the np.fft.ifft() function is the same as the actual Data. But I have no idea why.

Thank you for you time

import numpy as np
import matplotlib.pyplot as plt

X=np.arange(-0.6,0.6,0.00001)
original_FT=[]
FT_known_result=[]
for x_val in X:
    original_FT.append(1/(1 + x_val**2))
    FT_known_result.append(np.exp(-abs(x_val)))


FT_test_list=np.fft.ifftshift(original_FT)



plt.figure()
plt.plot(X,FT_test_list,label='FT calc ifft')
plt.plot(X,FT_known_result,label='real FT')
plt.plot(X,original_FT,label="orginal data")
plt.legend()
plt.show()

Solution

  • All of the below done in a Jupyter notebook. I didn't redo the analytical FT of a Lorentzian, so there may be an error in there. I'm assuming the following to be right:

    $$ FT  \left(  \frac{1}{\pi \Delta f} \frac{1} {1 + (\frac{f}{\Delta f})^2} \right)  = e^{ 2 \pi \Delta f |t|} $$
    

    enter image description here

    import numpy as np
    import matplotlib.pyplot as p
    %matplotlib inline
    
    
    start,stop,num=-.5,.5,1000
    t=np.linspace(start,stop,num)    # stop-start= 1.0 sec, inverse is frequency resolution , i.e 1 Hz
    dt=(stop-start)/num        
    print(f'time interval:  {dt} secs')     # will be 1 millisecond, inverse is frequency span i.e -500 Hz to +500 Hz
    
    freq = np.fft.fftfreq(num,dt)
    # print(freq)   #  0..499, -500.. -1  # unshifted; that can be useful, you don't need to shift to get the plots right
    freq2 =  np.fft.fftshift(freq)
    # print(freq2)    #  -500.. 0 .. 499    # we need the shift because you want to compare a spectrum = f (frequency)         
    print(f'freq interval :  {freq2[1] - freq2[0]}  Hz')                 
    # in numpy , given the time and frequency vectors we just made it is easy to get a function filled w/o a for loop
    # let's assume that original_FT is a spectrum as is typical for a Lorentzian
    
    df=5  # width of lorentzian
    
    ori= 1/np.pi /df /(1 + (freq2/(df ))**2)   # a lorentzian of variable width, FT of bisided exponential
    print(f'integral over lorentzian {np.sum(ori):.3f}') 
    analytical= np.exp(-abs(t*df*2*np.pi ))     # a bisided exponential  (the envelope of a time dependency), 
                                        # the width has to be inversely proportional to the width of the Lorentzian
    computed= np.abs(np.fft.fftshift(np.fft.fft(ori))) 
    
    
    p.figure(figsize=(10,4))
    p.subplot(131)
    p.plot(freq2,ori,label="orginal spectral data")
    p.xlabel('frequencies/Hz')
    p.title('freq. domain')
    p.legend()
    
    p.subplot(132)
    p.plot(t ,analytical ,label='real FT')
    p.plot(t ,computed ,label='calc ifft')
    p.xlabel('time/secs')
    p.title('time domain')
    p.legend()
    
    p.subplot(133)
    p.plot(t[490:510],analytical[490:510],'.-',label='real FT')
    p.plot(t[490:510],computed[490:510],'.-',label='calc ifft')
    p.xlabel('time/secs')
    p.title('zoomed in')
    p.legend();
    

    enter image description here