numpy.fft.fft() implementation in Python

Question

I am trying to obtain the spectrum of periodic signals using the fft function. Then plot the magnitude and phase of the transform. Magnitude plots are OK, but phase plots are completely unexpected. For example I used the functions Sin³(t) and Cos³(t). The code I have used is:

import matplotlib.pyplot as plt
import numpy.fft as nf
import math
import numpy as np
pi = math.pi

N=512

# Sin³(x)
t=np.linspace(-4*pi,4*pi,N+1);t=t[:-1]
y=(np.sin(t))**3
Y=nf.fftshift(nf.fft(y))/N
w=np.linspace(-64,64,N+1);w=w[:-1]

plt.figure("0")
plt.subplot(2,1,1)
plt.plot(w,abs(Y),'ro',lw=2)
plt.xlim((-4,4))
plt.ylabel(r"$|Y|$",size=16)
plt.title("Spectrum of sin\u00B3(t)")
plt.grid(True)
plt.subplot(2,1,2)
ii=np.where(abs(Y)>1e-3)
plt.plot(w[ii],np.angle(Y[ii]),'go',lw=2)
plt.xlim((-4,4))
plt.ylabel(r"Phase of $Y$",size=16)
plt.xlabel(r"$\omega$",size=16)
plt.grid(True)

# Cos³(x)

t=np.linspace(-4*pi,4*pi,N+1);t=t[:-1]
y=(np.cos(t))**3
Y=nf.fftshift(nf.fft(y))/N
w=np.linspace(-64,64,N+1);w=w[:-1]

plt.figure("1")
plt.subplot(2,1,1)
plt.plot(w,abs(Y),'ro',lw=2)
plt.xlim((-4,4))
plt.ylabel(r"$|Y|$",size=16)
plt.title("Spectrum of cos\u00B3(t)")
plt.grid(True)
plt.subplot(2,1,2)
ii=np.where(abs(Y)>1e-3)
plt.plot(w[ii],(np.angle(Y[ii])),'go',lw=2)
plt.xlim((-4,4))
plt.ylabel(r"Phase of $Y$",size=16)
plt.xlabel(r"$\omega$",size=16)
plt.grid(True) 

The plots obtained are:

i) For Sin³(t) - Magnitude and Phase of Spectrum of Sin³(t)

ii) For Cos³(t) - Magnitude and Phase of Spectrum of Cos³(t)

As You can see in the links above, The magnitude of Both functions are fine. The phase of Spectrum of Sin³(t) is correct as expected.

Since Cos³(t) is real and even, and expanding in terms of complex exponentials reveal that the Phase of the Coefficients is 0. But the plotted graph shows entirely different answer (see 2): at w=-3 the phase is near 5 radians. What is the mistake I have done and what is the correct way to implement the FFT.

Answer 1

np.fft is behaving as expected; it's your plot which is causing the confusion. A call to plt.tight_layout() should help clear things up for you:

If you take a careful look at the y-axis of the phase for cos³(t) you'll see that the values all have a 1e-17 prefactor. So at w = -3 the phase is not near 5 radians, it's actually near 5e-17 radians (which is for all intents and purposes zero).

I tidied up your code whilst I was looking into this:

import matplotlib.pyplot as plt
import numpy.fft as nf
import numpy as np

plt.rcParams['axes.grid'] = True

pi = np.pi
N = 512

t = np.linspace(-4*pi, 4*pi, N+1)[:-1]


def fft_func(x, func, N):
    y = func(x) ** 3
    Y = nf.fftshift(nf.fft(y)) / N
    w = np.linspace(-64, 64, N+1)[:-1]
    return y, Y, w


def plot_mag_phase(w, Y, func_name):
    fig, ax = plt.subplots(2, 1)

    ii = np.where(abs(Y) > 1e-3)

    ax[0].plot(w, abs(Y), 'ro', lw=2)
    ax[1].plot(w[ii], np.angle(Y[ii]), 'go', lw=2)

    ax[0].set_xlim(-4, 4)
    ax[1].set_xlim(-4, 4)

    ax[0].set_ylabel("$|Y|$", size=16)
    ax[1].set_ylabel("Phase of $Y$", size=16)

    ax[0].set_title("Spectrum of {}\u00B3(t)".format(func_name))
    ax[1].set_xlabel(r'$\omega$', size=16)

    fig.tight_layout()

# sin x ^ 3
y, Y, w = fft_func(t, np.sin, N)
plot_mag_phase(w, Y, 'sin')

# cos x ^ 3
y, Y, w = fft_func(t, np.cos, N)
plot_mag_phase(w, Y, 'cos')