Why numpy fft return incorrect phase information?

Question

I compare phase and amplitude spectrum in Matlab and numpy. I think Matlab work correct, but numpy compute correct amplitude spectrum, but phase spectrum is strange. How i must change python code for correct computing fft by numpy?

Matlab:

fs = 1e4;
dt = 1 / fs;
t = 0:dt:0.5;
F = 1e3;

y = cos(2*pi*F*t);
S = fftshift(fft(y) / length(y));
f_scale = linspace(-1, 1, length(y)) * (fs / 2);

a = abs(S);
phi = (angle(S));

subplot(2, 1, 1)
plot(f_scale, a)
title('amplitude')

subplot(2, 1, 2)
plot(f_scale, phi)
title('phase')

Python:

import numpy as np
import matplotlib.pyplot as plt

fs = 1e4
dt = 1 / fs
t = np.arange(0, 0.5, dt)
F = 1e3

y = np.cos(2*np.pi*F*t)
S = np.fft.fftshift(np.fft.fft(y) / y.shape[0])
f_scale = np.linspace(-1, 1, y.shape[0]) * (fs / 2)

a = np.abs(S)
phi = np.angle(S)

plt.subplot(2, 1, 1, title="amplitude")
plt.plot(f_scale, a)

plt.subplot(2, 1, 2, title="phase")
plt.plot(f_scale, phi)

plt.show()

matlab output

numpy output

Answer 1

It's a problem in understanding np.arange. It stops one dt before reaching the desired value (the interval you pass is open on the right side). If you define

t = np.arange(0, 0.5+dt, dt)

everything will work fine.

Answer 2

As pointed out in another answer, to make the Python plot match the matlab output, you have to adjust the t array to have the same values as the t array in the matlab code.

However, if your intent was to have an integer number of periods in the signal, so the FFT has just two nonzero values (at ± the input frequency), then it is the Python code that is correct. The phase in the Python code looks strange because all the Fourier coefficients except those associated with the signal's frequency are (theoretically) 0. With finite precision arithmetic, the coefficients end up being numerical "noise" with very small amplitude and essentially random phase.