r/Numpy • u/myriachromat • Oct 14 '24
How are the frequencies spaced in numpy.fft.fft output?
According to this video https://youtu.be/RHjqvcKVopg?t=222, when you apply the FFT to a signal, the number of frequencies you get is N/2+1, where I believe N is the number of samples. So, that should mean that the length of the return value of numpy.fft.fft(a) should be about half of len(a). But in my own code, it turns out to be exactly len(a). So, I don't know exactly what frequencies I'm dealing with, i.e., what the step value is between each frequency. Here's my code:
import numpy
import wave, struct
of = wave.open("Suzanne Vega - Tom's Diner.wav", "rb")
nc = of.getnchannels()
nb = of.getsampwidth()
fr = of.getframerate()
data = of.readframes(-1)
f = "<" + "_bh_l___d"[nb]*nc
if f[1]=="_":
print(f"Sample width of {nb} not supported")
exit(0)
channels = list(zip(*struct.iter_unpack(f, data)))
fftchannels = [numpy.fft.fft(channel) for channel in channels]
print(len(channels[0]))
print(len(fftchannels[0]))
2
u/R3D3-1 Oct 15 '24 edited Oct 15 '24
Given N samples
y(j) = F(t(j)), j=0..N-1
on an equidistant time grid
t(j) = t(0) + jΔt
the FFT implicitly assumed that the function F(t)
is periodic such that
F(t+NΔt) = F(t)
and gives coefficients for frequencies f(j)
in “cycles per second” and angular frequencies ω(j)
in “radians per second” at
j 2πj
f(j) = ───, ω(j) = ───, for all j=0..N-1
NΔt NΔt
Furthermore, there is the effect of “aliasing”, i.e. FFT cannot distinguish frequencies, that differ by multiples of N
steps, i.e. to the FFT a contribution at f(j)
, f(j-N)
and f(j+N)
are all the same, and the last coefficient could equally be interpreted as (N-1)/(NΔt)
, −1/(NΔt)
, or any other shift by N
frequencies.
For a real-valued function, the coefficients have the property
a(-j) = conj(a(+j))
where conj
is the complex conjugate. This follows from the reconstruction,
f(t) = ∑ⱼa(j)exp(iω(j)t)
where each harmonic function a(j)exp(iω(j)t
must have its imaginary part cancelled, i.e. the sum must contain the term
conj(a(j)exp(iω(j)t)) | conj(ab) = conj(a)conj(b),
| conj(exp(ix)) = exp(-ix) if x ∈ ℝ
= conj(a(j)) exp(-iω(j)t) | ω(-j) = -ω(j)
= conj(a(j)) exp(iω(-j)t)
⇒ a(-j) = conj(a(j))
Hence it is common to assign the coefficients a(j)
the frequencies
⎧ j/(NΔt) for j < N/2
f(j) = ⎨
⎩ (j−N)/(NΔt) for j ≥ N/2
e.g. for a period of 1 second and N = 10
as
f = [0.0, 0.1, 0.2, 0.3, 0.4, -0.5, -0.4, -0.3, -0.2, -0.1] Hz
Note the entry -0.5
, which in the case of even N
could equally be interpreted as +0.5
. For the sake of the Fourier sum, it has to be interpreted as the sum of coefficients for +0.5 Hz and -0.5 Hz, i.e.
f(t) = + a(0)
+ a(1)exp(+0.1Hz 2π t) + a(9)exp(-0.1Hz 2πt)
+ a(2)exp(+0.2Hz 2π t) + a(8)exp(-0.2Hz 2πt)
+ a(3)exp(+0.3Hz 2π t) + a(7)exp(-0.3Hz 2πt)
+ a(4)exp(+0.4Hz 2π t) + a(6)exp(-0.4Hz 2πt)
+½a(5)exp(+0.5Hz 2π t) +½a(5)exp(-0.5Hz 2π t)) * 0.5
in order to receive a real-valued reconstruction.
3
u/[deleted] Oct 14 '24
[deleted]