Output of fft.fft() for magnitude and phase (angle) not corresponding the the values set up

I set up a sine wave of a certain amplitude, frequency and phase, and tried recovering the amplitude and phase:

import numpy as np
import matplotlib.pyplot as plt

N = 1000  # Sample points     
T = 1 / 800         # Spacing

t = np.linspace(0.0, N*T, N) # Time
frequency = np.fft.fftfreq(t.size, d=T) # Normalized Fourier frequencies in spectrum.

f0 = 25              # Frequency of the sampled wave
phi = np.pi/6       # Phase
A = 50              # Amplitude
s = A * np.sin(2 * np.pi * f0 * t - phi) # Signal

S = np.fft.fft(s)   # Unnormalized FFT

fig, [ax1,ax2] = plt.subplots(nrows=2, ncols=1, figsize=(10, 5))
ax1.plot(t,s,'.-', label='time signal')  
ax2.plot(freq[0:N//2], 2/N * np.abs(S[0:N//2]), '.', label='amplitude spectrum')

plt.show()

index, = np.where(np.isclose(frequency, f0, atol=1/(T*N))) # Getting the normalized frequency close to f0 in Hz)
magnitude = np.abs(S[index[0]]) # Magnitude
phase = np.angle(S[index[0]])   # Phase
print(magnitude)
print(phase)
phi
#21785.02149316858
#-1.2093259641890741
#0.5235987755982988

Now the amplitude should be 50, instead of 21785, and the phase pi/6=0.524, instead of -1.2.

Am I misinterpreting the output, or the answer on the post referred to in the link above?

• You need to normalize the fft by 1/N with one of the two following changes (I used the 2nd one):
  S = np.fft.fft(s) --> S = 1/N*np.fft.fft(s)
magnitude = np.abs(S[index[0]]) --> magnitude = 1/N*np.abs(S[index[0]])
• Don't use index, = np.where(np.isclose(frequency, f0, atol=1/(T*N))), the fft is not exact and the highest magnitude may not be at f0, use np.argmax(np.abs(S)) instead which will give you the peak of the signal which will be very close to f0
• np.angle messes up (I think its one of those pi,pi/2 arctan offset things) just do it manually with np.arctan(np.real(x)/np.imag(x))
• use more points (I made N higher) and make T smaller for higher accuracy
• since a DFT (discrete fourier transform) is double sided and has peak signals in both the negative and positive frequencies, the peak in the positive side will only be half the actual magnitude. For an fft you need to multiply every frequency by two except for f=0 to acount for this. I multiplied by 2 in magnitude = np.abs(S[index])*2/N

N = 10000
T = 1/5000
...
index = np.argmax(np.abs(S))
magnitude = np.abs(S[index])*2/N
freq_max = frequency[index]
phase = np.arctan(np.imag(S[index])/np.real(S[index])) 
print(f"magnitude: {magnitude}, freq_max: {freq_max}, phase: {phase}") print(phi)

Output: magnitude: 49.996693276663564, freq_max: 25.0, phase: 0.5079341239733628