Understanding the FFT phase spectrum with a simple example

Question

I'm trying to compute the DFT using scipy's functions. I don't understand why the phase spectrum of a simple sine wave with 2 Hz frequency doesn't show $\pm\pi/2$ at the $\pm 2Hz$ frequencies. Instead, the phase plot seems to have some linear dependence in the frequency, which I don't understand. I provide the code for assistance. How can this be fixed? Looks like a simple issue I'm not grasping. Please help.

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft, fftfreq

# 1Hz sine wave
npts  = 100
tmax  = 10                          
t     = np.linspace(0, tmax, npts)
y     = np.sin(2*np.pi*2*t)
dt    = tmax/npts

# FFT computation
Y      = fft(y)
freq   = fftfreq(npts, d=dt)
amplit = abs(Y)/npts
phase  = np.angle(Y)
phase  = phase / np.pi

fig, ax = plt.subplots(1, 2, figsize=(10,4))
ax[0].plot(freq, amplit)
ax[1].plot(freq, phase)

ax[0].plot([2,2], [0,0.5], '--')
ax[0].plot([-2,-2], [0,0.5], '--')
ax[1].plot([2,2], [-1,1], '--')
ax[1].plot([-2,-2], [-1,1], '--')

Answer 1

For sinusoids that are not exactly integer periodic in the FFT length, an FFT measures the phase at a circular discontinuity. And that discontinuity flips direction as frequency changes from slightly below to slightly above an exact integer periodic-in-aperture frequency. This is part of the effect of the default rectangular windowing of any finite length DFT or FFT. Thus making the FFT phase result of a frequency sweep hard to interpret.

The way to remove this measurement at a circular discontinuity is to shift the waveform discontinuity away from the phase measurement point (FFT input sample 0) by doing a (circular) FFTShift before the FFT. Then your FFT phase result will be referenced to the center of the original input vector (original waveform sample N/2), which will be continuous for continuous input functions or waveforms.

After the FFT following an fftshift, you can use the FFT frequency and phase results for the center of the original input waveform to compute the phase at any other point (for instance the beginning, or sample 0 of the original input waveform. Or the end, etc.).

Answer 2

As @hotpaw2 says, the FFT assumes periodic boundary conditions, that is the next point on the right should be equal to the first point on the left. It's not the case in your code because the time $t=10$ is included in your time array. The next point on the right is thus $t=10+dt$ and you get the discontinuity because $y(10+dt) \neq y(0)$. You can fix this by defining the time array as t=np.linspace(0, tmax, npts, endpoint=False). Now this fix will only work as long as the time span is an integer number of the signal period. If it's not you'll have to rely on @hotpaw2's more general solution.

PS: if your version of scipy allows it, you should use the newer submodule scipy.fft instead of the legacy submodule scipy.fftpack.