Discrete Fourier Transform in Python

[Silviu]

Homework Statement

I need to calculate the derivative of a function using discrete Fourier transform (DFT). Below is a simplified version of my code (just for sin function) in python

Homework Equations

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)

The Attempt at a Solution

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!

 

[collinsmark]

Homework Statement

I need to calculate the derivative of a function using discrete Fourier transform (DFT). Below is a simplified version of my code (just for sin function) in python

Homework Equations

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)

The Attempt at a Solution

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!

========== For lurkers reading this thread, here is some background on what is trying to be done. ===
https://en.wikibooks.org/wiki/Paral...ng_Derivatives_using_Fourier_Spectral_Methods
=========================================================================

The problem lies in how you are choosing your "k" values. As you have them above, they are increasing linearly on the segment [0, N), where N is the number of samples. In other words, in this case since N = 4, your "k" values are [0, 1, 2, 3]. But that's not what you want.

Instead, what you want is for your "k" values to increase linearly along the segment [0, N/2), and then the second half be from [-N/2, 0).

Let me explain by example. In your present implementation you have 4 sample points. So you want your array of "k" values to be [0, 1, -2, -1]

If you instead were using 8 sample points, you want you "k" values to be [0, 1, 2, 3, -4, -3, -2, -1]

-----

In your program, if you want to generalize this to use more sample points than just 4, define "N" somewhere to be the number of sample points you want. Then define your "theta" and "k" parameters to be a function of N (rather than hardcoding the number "4" as you did above).

You can properly make your array of "k" values by utilizing Python's numpy method numpy.concatenate(). Make the upper part of the array and the lower part of the array separately, both using np.arange(). Then you can combine them into a single array using np.concatenate().

Something like (with details left out)
a = np.arange( ''' lower part of k array ''')
b = np.arange( ''' upper part of k array ''')
k = np.concatenate((a, b), axis = 0)

[Edit: the above advice assumes that the number of samples, N, is even. If N is odd, there is the special case for the center element of the "k" array that you need to handle separately.]

 

[Silviu]

========== For lurkers reading this thread, here is some background on what is trying to be done. ===
https://en.wikibooks.org/wiki/Paral...ng_Derivatives_using_Fourier_Spectral_Methods
=========================================================================

The problem lies in how you are choosing your "k" values. As you have them above, they are increasing linearly on the segment [0, N), where N is the number of samples. In other words, in this case since N = 4, your "k" values are [0, 1, 2, 3]. But that's not what you want.

Instead, what you want is for your "k" values to increase linearly along the segment [0, N/2), and then the second half be from [-N/2, 0).

Let me explain by example. In your present implementation you have 4 sample points. So you want your array of "k" values to be [0, 1, -2, -1]

If you instead were using 8 sample points, you want you "k" values to be [0, 1, 2, 3, -4, -3, -2, -1]

-----

In your program, if you want to generalize this to use more sample points than just 4, define "N" somewhere to be the number of sample points you want. Then define your "theta" and "k" parameters to be a function of N (rather than hardcoding the number "4" as you did above).

You can properly make your array of "k" values by utilizing Python's numpy method numpy.concatenate(). Make the upper part of the array and the lower part of the array separately, both using np.arange(). Then you can combine them into a single array using np.concatenate().

Something like (with details left out)
a = np.arange( ''' lower part of k array ''')
b = np.arange( ''' upper part of k array ''')
k = np.concatenate((a, b), axis = 0)

[Edit: the above advice assumes that the number of samples, N, is even. If N is odd, there is the special case for the center element of the "k" array that you need to handle separately.]

Thank you for your reply. However the values of k are given in the problem. Even the plot of $c_k$ vs k is there so I can't really change that...

 

[phyzguy]

Thank you for your reply. However the values of k are given in the problem. Even the plot of $c_k$ vs k is there so I can't really change that...

collinsmark is right. You must be misinterpreting the problem. Can you share the exact problem statement you were given?

 

[Silviu]

collinsmark is right. You must be misinterpreting the problem. Can you share the exact problem statement you were given?

I attached it. I apologize if I misinterpret something.

 

[phyzguy]

The point is that both the positive and negative frequencies need to be included. For real inputs like you have, the positive and negative frequencies contain the same information, so you can ignore the negative frequencies and just plot ck for the positive frequencies. But you need to include the negative frequencies in the k array when you use it to calculate the derivative. You might try reading this description of the numpy.fft function. I've quoted the key passage below. Note that there is the np.fft.fftfreq() function that will return the k values.

"The values in the result(they mean here the k values) follow so-called “standard” order: If A = fft(a, n), then A[0] contains the zero-frequency term (the sum 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. For an even number of input points, A[n/2] represents both positive and negative Nyquist frequency, and is also purely real for real input. For an odd number of input points, A[(n-1)/2] contains the largest positive frequency, while A[(n+1)/2] contains the largest negative frequency. The routine np.fft.fftfreq(n) returns an array giving the frequencies of corresponding elements in the output. The routine np.fft.fftshift(A) shifts transforms and their frequencies to put the zero-frequency components in the middle, and np.fft.ifftshift(A) undoes that shift."

 

[Silviu]

The point is that both the positive and negative frequencies need to be included. For real inputs like you have, the positive and negative frequencies contain the same information, so you can ignore the negative frequencies and just plot ck for the positive frequencies. But you need to include the negative frequencies in the k array when you use it to calculate the derivative. You might try reading this description of the numpy.fft function. I've quoted the key passage below. Note that there is the np.fft.fftfreq() function that will return the k values.

"The values in the result(they mean here the k values) follow so-called “standard” order: If A = fft(a, n), then A[0] contains the zero-frequency term (the sum 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. For an even number of input points, A[n/2] represents both positive and negative Nyquist frequency, and is also purely real for real input. For an odd number of input points, A[(n-1)/2] contains the largest positive frequency, while A[(n+1)/2] contains the largest negative frequency. The routine np.fft.fftfreq(n) returns an array giving the frequencies of corresponding elements in the output. The routine np.fft.fftshift(A) shifts transforms and their frequencies to put the zero-frequency components in the middle, and np.fft.ifftshift(A) undoes that shift."

Ah, ok I understand what you mean now. However, for the first part of my question, why does that decay point depends on n, I don't need to negative frequencies. So, is my explanation correct for that? Thank you!

 

[phyzguy]

Ah, ok I understand what you mean now. However, for the first part of my question, why does that decay point depends on n, I don't need to negative frequencies. So, is my explanation correct for that? Thank you!

Which explanation?