2

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?

Hersh Joshi
  • 419
  • 3
  • 13
Antoni Parellada
  • 4,253
  • 6
  • 49
  • 114
  • 1
    Your phase is correct. A cosine without shift has a phase of 0. A sine without a shift will have a phase of -pi/2. Your signal should have a phase of -pi/2+pi/6. – Cris Luengo Jan 19 '22 at 00:56

1 Answers1

3
  • 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

Antoni Parellada
  • 4,253
  • 6
  • 49
  • 114
Hersh Joshi
  • 419
  • 3
  • 13
  • THANK YOU!!! (upper case intended!) I just read the answer quickly but I see the output being as expected! There is justice in the world :-) (after I read carefully if I don't have any questions I'll come back to accept). Thanks again! – Antoni Parellada Jan 18 '22 at 17:48
  • Do you know why if I change the phase to `np.pi/2` the phase extracted is numerically zero - I understand I'm basically turning the cosine into a sine... – Antoni Parellada Jan 18 '22 at 18:04
  • Likewise, if I change the phase to phase: `phi = np.pi/8 ~ 0.39 ` I get `magnitude: 50.000000000000014, freq_max: -25.0, phase: 1.178097245096179`. – Antoni Parellada Jan 18 '22 at 18:07
  • `frequency` is both positive and negative. The absolute value of an fft will have peaks in both the positive and negative. These peaks may differ in phase, but their absolute value is mathematically identical (although the numerical calculation may be slightly different) and you don't need to worry about it. for `phi=np.pi/2` it gives me `1.55` as the output so I'm unsure what you mean – Hersh Joshi Jan 18 '22 at 18:14
  • if you look at the DFT formula, https://en.wikipedia.org/wiki/Discrete_Fourier_transform. k is the index and is evaluated on [-N/2,N/2]. Giving negative frequencies (its also sometimes evaluated on [0,N] instead but it doesnt look like numpy is doing that). Also see: https://stackoverflow.com/questions/49172828/wrong-amplitude-of-the-fft-of-a-signal – Hersh Joshi Jan 18 '22 at 18:22
  • Here is the output I get with the code after your answer for `phi = np.pi/2`: `magnitude: 50.0, freq_max: 25.0, phase: -1.2582089486833056e-14`. So the phase is practically zero, instead of 1.57. If you get around the correct value when you run it, could you do me a big favor and just include the code for that value? – Antoni Parellada Jan 18 '22 at 18:28
  • Thank you for touching on it, but I do understand the mirror image of the FT for real signal with half the magnitude on either side, and the need for normalization, which you implemented beautifully. Not so sure about the relationship to the phase offset, though. – Antoni Parellada Jan 18 '22 at 18:30
  • I understand what the problem is: I think you made a mistake transcribing the code in your answer. The arctan should be imaginary/real, not the way you have it. Also there is a parenthesis issue. It should be like follows: `phase = np.arctan(np.imag(S[index])/np.real(S[index])) print(f"magnitude: {magnitude}, freq_max: {freq_max}, phase: {phase}") print(phi)` – Antoni Parellada Jan 18 '22 at 20:57
  • `arctan` gives you only the phase modulo pi. You need to use `atan2(y,x)` to get a proper phase. Also, in OP’s situation, the obtained phase was correct, 0.5 is not. – Cris Luengo Jan 19 '22 at 00:53