1

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()
Scott Stensland
  • 26,870
  • 12
  • 93
  • 104
Stefan
  • 23
  • 3
  • You are not using the FFT, you are just using the shift associated with the FFT. Just add `y = np.fft.fft(original_FT)`. One should also sample the function to be transformed well into the wings (i.e. where the values become small) or one will have significant aliasing (i.e. make your function narrower if you want to keep the range). Also the number of datapoints is not needed to get the picture. Lastly, you can't plot the `time` and `frequency` results in the same plot. They do have the same number of points, but different physical meaning. – roadrunner66 Jan 09 '20 at 03:05
  • Thank you for the quick response. I still have a question though: Do you mean I have to add y= np.fft.fft(original FT) and then invert it again? Since I am not looking for the FT but the inverse FT. I have found an error since the vector must have a certain structure. But again the Data does hase a certain variance to the real FT, which I find odd. Regarding the time and freq. yes you are 100 correct that those are not the same I just added it into the plot for visualization reasons. – Stefan Jan 09 '20 at 10:12

1 Answers1

0

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

roadrunner66
  • 7,772
  • 4
  • 32
  • 38
  • Thank you for your answer. The key line that I was missing was np.fft.fftshift(np.fft.fft(ori)) – Stefan Jan 10 '20 at 13:30