Computational Science Asked by user37540 on June 10, 2021
I write the following fast Fourier transform code into my Python notebook expecting to see a plot wherein there’s a spike at $1/2pi$ since that’s the frequency of the sin function, but instead I get the plot below. Any ideas as to what’s going on?
import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.fft import fft
x = np.arange(1000)
y = np.sin(x)
freq = fft(y)
plt.plot(np.real(freq))
plt.show()
The base FFT is defined for both negative and positive frequencies. What you see here is not what you think. Scipy returns the bin of the FFT in that order: positive frequencies from 0 to fs/2, then negative frequencies from -fs/2 up to 0. If your signal is real, you should look at the real FFT:https://docs.scipy.org/doc/scipy/reference/generated/scipy.fft.rfft.html#scipy.fft.rfft .
Scipy's fft does only give you the bin height, not their frequencies. To obtain them, use fftfreq. To properly reorder the FFT bins, use scipy's fftshift: https://docs.scipy.org/doc/scipy/reference/generated/scipy.fftpack.fftshift.html#scipy.fftpack.fftshift
Note that numpy also has a similar FFT module.
Here is a quick adaptation of your code:
import math
import matplotlib.pyplot as plt
import numpy as np
import scipy.fft as fft
f=0.1 # input signal frequency Hz
T = 10*1/f # duration of the signal
fs = f*4 # sampling frequency (at least 2*f)
x = np.arange(0,T,1/fs) # time vector of the sampling
y = np.sin(2*np.pi*f*x) # sampled values
# compute the FFT bins, diving by the number of samples to get
# the true "physical" amplitudes
fft_values = fft.fft(y)/x.size
# determine the corresponding physical frequencies
freq = fft.fftfreq(x.size, x[1]-x[0])
# change the ordering so that frequencies are increasing
fft_values = fft.fftshift(fft_values)
freq = fft.fftshift(freq)
# plot the result
fig, ax = plt.subplots(2,1,sharex=True)
ax[0].plot(freq, np.abs(fft_values))
ax[0].set_ylabel('norm')
# plot the phase, only for the sufficiently large bins
ax[1].plot(freq, np.angle(fft_values)*( np.abs(fft_values) > 1e-2))
ax[1].set_ylabel('phase (rad)')
ax[-1].set_xlabel('f (Hz)')
ax[0].grid(); ax[1].grid()
And here is the result:
I would suggest reading the explanation here: https://docs.scipy.org/doc/scipy/reference/generated/scipy.fft.fft.html#scipy.fft.fft and https://docs.scipy.org/doc/scipy/reference/tutorial/fft.html#d-discrete-fourier-transforms
Why do you obtain two bins ? Your sinus input can be expressed as follows: $sin(2pi f t) = frac{exp(2 i pi f t) - exp(-2 i pi f t)}{2 i}$, therefore you have one negative and one positive frequencies in your FFT results, both of which are given only half of your sinusoid's amplitude. Their FFT amplitudes are actually pure imaginary, therefore the corresponding phases are $pm frac{pi}{2}$.
In your initial test, you were only considering the real part of the FFT. Had you considered a number of samples that is a multiple of the number of period in your signal, you would have obtained 0 ! It's only because you sampling frequency was not a multiple of your signal's frequency that you were able to have a (misleading) result because of a phenomenon called "spectral leakage".
(I haven't worked with an FFT in a while, so please feel free to correct my mistakes if you find some)
Answered by Laurent90 on June 10, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP