Discrete Fourier Transform in Python

Question

I need to use discrete Fourier transform (DFT) in Python (and inverse DFT) and the results I obtain are a bit weird, so I tried on a small example and I am not sure I understand the mistake (if it is math or coding). Here is my small version of the code:

from __future__ import division
import numpy as np
from pylab import *

pi = np.pi

def f(x):
    return sin(x)

theta = np.arange(0,2*pi,2*pi/4)
k = np.arange(0,4,1)

x = f(theta)
y = np.fft.fft(x)

derivative = np.fft.ifft(1j*k*y) 
print(derivative)

So what I do is to sample sin at 4 different points between 0 and 2pi and create with these numbers a vector x. Then I take the DFT of x to get y. What I want is to get the derivative of sin at the chosen points, so to do this I multiply y by k (the wave number, which in this case would be 0,1,2,3) and my the imaginary number 1j (this is because in the Fourier sum I have for each term something of the form e^{ikx}). So in the end I take the inverse DFT of 1jky and I am supposed to get the derivative of sin. But what I get is this.

[ -1.00000000e+00 -6.12323400e-17j  -6.12323400e-17 +2.00000000e+00j
   1.00000000e+00 +1.83697020e-16j   6.12323400e-17 -2.00000000e+00j]

when I was supposed to get this

[1,0,-1,0]

ignoring round-off errors. Can someone tell me what am I doing wrong? Thank you!

score 0 · Answer 1 · answered Oct 29 '17 at 01:26

Manipulation of the spectrum must preserve this Hermitian symmetry if the inverse FFT is to yield result. Accordingly, the derivative operator in the frequency domain is defined over the lower half of the spectrum, and the upper half of the spectrum constructed from symmetry. Note that for spectrum of even sizes the value at exactly N/2 must be its own symmetry, hence must have a imaginary part which is 0. The following illustrate how to construct this derivative operator:

N = len(y)
if N%2:
  derivative_operator = np.concatenate((np.arange(0,N/2,1),[0],np.arange(-N/2+1,0,1)))*1j
else:
  derivative_operator = np.concatenate((np.arange(0,N/2,1),np.arange(-N//2+1,0,1)))*1j

You'd use this derivative_operator in the frequency-domain as follow:

derivative = np.fft.ifft(derivative_operator*y)

In your sample case you should then get the following result

[  1.00000000e+00+0.j   6.12323400e-17+0.j  
  -1.00000000e+00+0.j  -6.12323400e-17+0.j]

which is within roundoff errors of your expected [1,0,-1,0].

Discrete Fourier Transform in Python

1 Answers1