20

Consider a function f(t), how do I compute the continuous Fouriertransform g(w) and plot it (using numpy and matplotlib)?

This or the inverse problem (g(w) given, plot of f(t) unknown) occurs if there exists no analytical solution to the Fourier Integral.

thomasfermi
  • 1,183
  • 1
  • 10
  • 20
  • +1 So you write a question and then answer it by yourself? – GingerHead Jun 06 '14 at 09:04
  • 13
    Yes, I read that people are encouraged to do so. This was one of the few numpy/matplotlib problems I didn't find a solution for by using google. So i thought I'd share the solution. The page where I read about answering your own question was here http://blog.stackoverflow.com/2011/07/its-ok-to-ask-and-answer-your-own-questions/ – thomasfermi Jun 06 '14 at 09:07
  • Hi, would you be so kind to have a look at my question [here](http://stackoverflow.com/questions/34428886/discrete-fourier-transfromation-from-a-list-of-x-y-points)? – Michael Kim Dec 23 '15 at 05:09

1 Answers1

32

You can use the numpy FFT module for that, but have to do some extra work. First let's look at the Fourier integral and discretize it: Here k,m are integers and N the number of data points for f(t). Using this discretization we get enter image description here

The sum in the last expression is exactly the Discrete Fourier Transformation (DFT) numpy uses (see section "Implementation details" of the numpy FFT module). With this knowledge we can write the following python script

import numpy as np
import matplotlib.pyplot as pl

#Consider function f(t)=1/(t^2+1)
#We want to compute the Fourier transform g(w)

#Discretize time t
t0=-100.
dt=0.001
t=np.arange(t0,-t0,dt)
#Define function
f=1./(t**2+1.)

#Compute Fourier transform by numpy's FFT function
g=np.fft.fft(f)
#frequency normalization factor is 2*np.pi/dt
w = np.fft.fftfreq(f.size)*2*np.pi/dt


#In order to get a discretisation of the continuous Fourier transform
#we need to multiply g by a phase factor
g*=dt*np.exp(-complex(0,1)*w*t0)/(np.sqrt(2*np.pi))

#Plot Result
pl.scatter(w,g,color="r")
#For comparison we plot the analytical solution
pl.plot(w,np.exp(-np.abs(w))*np.sqrt(np.pi/2),color="g")

pl.gca().set_xlim(-10,10)
pl.show()
pl.close()

The resulting plot shows that the script works enter image description here

thomasfermi
  • 1,183
  • 1
  • 10
  • 20
  • 4
    +1 For nice way of publishing a question and answering it by yourself! – GingerHead Jun 06 '14 at 09:03
  • If you're computing a (discreticized) FT this way, can you include frequencies smaller than 1 or periods that are multiples of the length of the input array? Perhaps you could help answering this question: http://dsp.stackexchange.com/questions/39017/looking-for-cycles-of-periods-longer-than-the-input-signal-length – mac13k Apr 09 '17 at 08:37
  • Yeah, for a frequency-to-time-Fourier-Transform you SHOULD include small frequencies, otherwise your result for long times will not be very good. I think your question is not directly related, and I cannot answer it without putting considerable research into it myself, sorry. – thomasfermi Apr 09 '17 at 10:22
  • I thought it was related, because someone suggested in the comments that the continuous transform could be the solution, because it can consider sines of periods longer than the input signal, while the discrete transform does not go beyond the signal length. I'm not if or how the discreticised continuous transform can do the same, but I will investigate... – mac13k Apr 09 '17 at 11:58
  • What is that green horizontal line? – rainman Jan 21 '19 at 15:54
  • @omehoque: That is because the frequency array w is not ordered. – thomasfermi Jan 28 '19 at 19:24
  • @thomasfermi: Thanks. I eventually figured it out. And thanks again for the excellent question and answer. – rainman Jan 28 '19 at 21:12
  • Not a Python question, but what if i do not want the time discretization symmetric about zero but with an offset, i.e. `t=offset + np.arange(t0,-t0,dt)`. What will the corresponding omegas be then? – Christoph90 May 22 '19 at 07:28
  • 1
    @Christoph90: If you do that, you will have the same frequencies and you will compute the Fourier transform of f(t-offset) instead of f(t) – thomasfermi Jun 12 '19 at 18:23