1

I have been thinking about it for a long time, but I don't find out what the problem is. Hope you can help me, Thank you.

F(s) Gaussian function
F(s)=1/(√2π s) e^(-(w-μ)^2/(2s^2 ))

Code:

import numpy as np
from matplotlib import pyplot as plt
from math import pi
from scipy.fft import fft

def F_S(w, mu, sig):
    return (np.exp(-np.power(w-mu, 2)/(2 * np.power(sig, 2))))/(np.power(2*pi, 0.5)*sig)

w=np.linspace(-5,5,100)
plt.plot(w, np.real(np.fft.fft(F_S(w, 0, 1))))
plt.show()

Result:

enter image description here

Dineshkarthik Raveendran
  • 726
  • 7
  • 19
user12807422
  • 23
  • 5
  • what do you expect it to be? – Marc Jan 29 '20 at 16:04
  • I want to get the curve after Fourier Transform of the Gaussian function – user12807422 Jan 29 '20 at 16:16
  • and how should that look compared to what you get as an answer? – Marc Jan 29 '20 at 16:17
  • 1
    Gaussian function should still be Gaussian after FT – user12807422 Jan 29 '20 at 16:18
  • If you want the plot to look like another Gaussian, you need to reorder coefficients, fix up the scalings, etc., according to whatever your DFT routine is doing. – bg2b Jan 29 '20 at 16:20
  • 1
    You should not take the real component of the result, but its absolute value – dedObed Jan 29 '20 at 16:22
  • Does this answer your question? [Analytical Fourier transform vs FFT of functions in Matlab](https://stackoverflow.com/questions/49317834/analytical-fourier-transform-vs-fft-of-functions-in-matlab) – Cris Luengo Jan 29 '20 at 17:21
  • The linked answer uses MATLAB, not Python, but the concepts are identical. It should explain exactly what the issue is with your code. – Cris Luengo Jan 29 '20 at 17:22
  • Thanks for everyone, the problem is solved, I found a related question https://stackoverflow.com/questions/5398304/fourier-transform-of-a-gaussian-is-not-a-gaussian-but-thats-wrong-python – user12807422 Jan 29 '20 at 17:40

4 Answers4

0

you have to change from time scale to frequency scale

Mohammed Khalid
  • 155
  • 1
  • 6
  • please look at these examples to see how to transform from time scale into frequency scale https://docs.scipy.org/doc/scipy-0.14.0/reference/tutorial/fftpack.html – Mohammed Khalid Jan 29 '20 at 17:35
0

When you make a FFT you will get the simetric tranformation, i.e, mirror of the positive to negative curve. Usually, you only will look at the positive side. Also, you should take care with sample rate, as FFT is designed to transform time domain input to frequency domain, the time, or sample rate, of input info matters. So add timestep in np.fft.fftfreq(n, d=timestep) for your sample rate.

If you simple want to make a fft of normal dist signal, here is another question with it and some good explanations on why are you geting this behavior:

Fourier transform of a Gaussian is not a Gaussian, but thats wrong! - Python

Guinther Kovalski
  • 1,629
  • 1
  • 7
  • 15
0

There are two mistakes in your code:

  • Don't take the real part, take the absoulte value when plotting.
  • From the docs:

If A = fft(a, n), then A[0] contains the zero-frequency term (the mean of the signal), which is always purely real for real inputs. Then A[1:n/2] contains the positive-frequency terms, and A[n/2+1:] contains the negative-frequency terms, in order of decreasingly negative frequency.

You can rearrange the elements with np.fft.fftshift.

The working code:

import numpy as np
from matplotlib import pyplot as plt
from math import pi
from scipy.fftpack import fft, fftshift

def F_S(w, mu, sig):
    return (np.exp(-np.power(w-mu, 2)/(2 * np.power(sig, 2))))/(np.power(2*pi, 0.5)*sig)

w=np.linspace(-5,5,100)
plt.plot(w, fftshift(np.abs(np.fft.fft(F_S(w, 0, 1)))))
plt.show()

Also, you might want to consider scaling the x axis too.

Péter Leéh
  • 2,069
  • 2
  • 10
  • 23
0

As was mentioned before you want the absolute value, not the real part. A minimal example, showing the the re/im , abs/phase spectra.

import numpy as np
import matplotlib.pyplot as p
%matplotlib inline

n=1001                       # add 1 to keep the interval a round number when using linspace
t = np.linspace(-5, 5, n )   # presumed to be time
dt=t[1]-t[0]                 # time resolution
print(f'sampling every  {dt:.3f} sec , so at {1/dt:.1f} Sa/sec, max. freq will be {1/2/dt:.1f} Hz')


y = np.exp(-(t**2)/0.01)      # signal in time

fr= np.fft.fftshift(np.fft.fftfreq(n, dt))  # shift helps with sorting the frequencies for better plotting
ft=np.fft.fftshift(np.fft.fft(y))           # fftshift only necessary for plotting in sequence

p.figure(figsize=(20,12))
p.subplot(231)
p.plot(t,y,'.-')
p.xlabel('time (secs)')
p.title('signal in time')

p.subplot(232)
p.plot(fr,np.abs(ft), '.-',lw=0.3)      
p.xlabel('freq (Hz)')
p.title('spectrum, abs');

p.subplot(233)
p.plot(fr,np.real(ft), '.-',lw=0.3)     
p.xlabel('freq (Hz)')
p.title('spectrum, real');

p.subplot(235)
p.plot(fr,np.angle(ft), '.-', lw=0.3)    
p.xlabel('freq (Hz)')
p.title('spectrum, phase');

p.subplot(236)
p.plot(fr,np.imag(ft), '.-',lw=0.3)      
p.xlabel('freq (Hz)')
p.title('spectrum, imag');

enter image description here

roadrunner66
  • 7,772
  • 4
  • 32
  • 38