5

I have a signal for which I need to calculate the magnitude and phase at 200 Hz frequency only. I would like to use Fourier transform for it. I am very new to signal processing. And this is my first time using a Fourier transform.

I found that I can use the scipy.fftpack.fft to calculate the FFT of the signal. Then use numpy.mag and numpyh.phase to calculate the magnitude and phases of the entire signal. But I would like to get the magnitude and phase value of the signal corresponding to 200 Hz frequency only. How can I do this using Python?

So far I have done.

from scipy.fftpack import fft
import numpy as np

fft_data = fft(signal)
magnitude = np.mag(fft_data)
phase = np.phase(fft_data)
Gilles
  • 3,386
  • 3
  • 21
  • 28
thileepan
  • 183
  • 1
  • 1
  • 8

1 Answers1

6

You can find the index of the desired (or the closest one) frequency in the array of resulting frequency bins using np.fft.fftfreq function, then use np.abs and np.angle functions to get the magnitude and phase.

Here is an example using fft.fft function from numpy library for a synthetic signal.

import numpy as np
import matplotlib.pyplot as plt

Number of sample points


N = 1000


Sample spacing


T = 1.0 / 800.0     # f = 800 Hz


Create a signal


x = np.linspace(0.0, NT, N)
t0 = np.pi/6   # non-zero phase of the second sine
y = np.sin(50.0  2.0np.pix) + 0.5np.sin(200.0  2.0np.pix + t0)
yf = np.fft.fft(y) # to normalize use norm='ortho' as an additional argument


Where is a 200 Hz frequency in the results?


freq = np.fft.fftfreq(x.size, d=T)
index, = np.where(np.isclose(freq, 200, atol=1/(T*N)))


Get magnitude and phase


magnitude = np.abs(yf[index[0]])
phase = np.angle(yf[index[0]])
print("Magnitude:", magnitude, ", phase:", phase)


Plot a spectrum


plt.plot(freq[0:N//2], 2/N*np.abs(yf[0:N//2]), label='amplitude spectrum')   # in a conventional form
plt.plot(freq[0:N//2], np.angle(yf[0:N//2]), label='phase spectrum')
plt.legend()
plt.grid()
plt.show()

And here is a useful manual with detailed explanations: reference.

megasplash
  • 396
  • 2
  • 13
  • thank you for the clear explanation. I have two doubts. 1. What should I do if I want to change the phase of the signals to 30 degrees? 2. Why are you not scaling the magnitude values with (2/N)? – thileepan Dec 18 '20 at 11:11
  • 1
    @thileepan 1. The phase t0 would be an additional term in the argument of a sine: A*sin(wt+t0). t0 = np.pi/6 should shift the signal to 30 degrees. 2. The example shows the default fft results. You can normalize the magnitude by setting the "norm" parameter like this: yf = np.fft.fft(y, norm='ortho'). Btw, my bad, np.isclose does not work as intended. I will fix it in the answer. – megasplash Dec 18 '20 at 14:03
  • @thileepan Regarding the normalization of the results, you may find this question to be helpful. – megasplash Dec 18 '20 at 14:33
  • Silly question (probably): I ran the code, and the output for the phase is phase: -0.2672456406675483, but you had set it up to be t0 = np.pi/6 # non-zero phase of the second sine. This is 0.5235987755982988. Why are they not the same? – Antoni Parellada Jan 18 '22 at 00:36
  • @AntoniParellada No stupid questions on this site) If you look at the phase spectrum plot, you can see that non-zero values are spread through all spectrum. That is why you see the difference. Also you can see that amplitude plot also has non-zero values through all spectrum. This happens when sine frequencies are not the exact multiple of the fundamental frequency of the DFT vector space, which was 2pi/1000 in my example. – megasplash Feb 01 '22 at 20:13
  • @megasplash Thank you! – Antoni Parellada Feb 01 '22 at 23:05
  • @megasplash, but why the phase spectrum spreads all the frequency range? why it is not just a point corresponding to the frequency of interest? – Curious May 16 '22 at 17:20
  • 1
    @Curious Calculating DFT means that you represent the signal as a linear combination of multiples of a fundamental frequency (0, f, 2f, 3f, 4f, ... , (N-1)f). When a sinusoid frequency in a signal is not an exact multiple of f it will get a contribution from every DFT coefficient to compensate that, which looks like non-zero values through the spectrum. This applies to both real and imaginary parts, i.e. to both amplitude and phase. – megasplash May 24 '22 at 11:19
  • @megasplash, and this explains such a phase incursion during the whole frequency range? – Curious May 24 '22 at 11:38
  • 1
    @Curious This explains the expected phase values spread. There is also a numerical computation error, which will alter the phase values even if an amplitude has single non-zero values. – megasplash May 24 '22 at 15:33