54

I used fft function in numpy which resulted in a complex array. How to get the exact frequency values?

3 Answers 3

75

np.fft.fftfreq tells you the frequencies associated with the coefficients:

import numpy as np

x = np.array([1,2,1,0,1,2,1,0])
w = np.fft.fft(x)
freqs = np.fft.fftfreq(len(x))

for coef,freq in zip(w,freqs):
    if coef:
        print('{c:>6} * exp(2 pi i t * {f})'.format(c=coef,f=freq))

# (8+0j) * exp(2 pi i t * 0.0)
#    -4j * exp(2 pi i t * 0.25)
#     4j * exp(2 pi i t * -0.25)

The OP asks how to find the frequency in Hertz. I believe the formula is frequency (Hz) = abs(fft_freq * frame_rate).

Here is some code that demonstrates that.

First, we make a wave file at 440 Hz:

import math
import wave
import struct

if __name__ == '__main__':
    # http://stackoverflow.com/questions/3637350/how-to-write-stereo-wav-files-in-python
    # http://www.sonicspot.com/guide/wavefiles.html
    freq = 440.0
    data_size = 40000
    fname = "test.wav"
    frate = 11025.0
    amp = 64000.0
    nchannels = 1
    sampwidth = 2
    framerate = int(frate)
    nframes = data_size
    comptype = "NONE"
    compname = "not compressed"
    data = [math.sin(2 * math.pi * freq * (x / frate))
            for x in range(data_size)]
    wav_file = wave.open(fname, 'w')
    wav_file.setparams(
        (nchannels, sampwidth, framerate, nframes, comptype, compname))
    for v in data:
        wav_file.writeframes(struct.pack('h', int(v * amp / 2)))
    wav_file.close()

This creates the file test.wav. Now we read in the data, FFT it, find the coefficient with maximum power, and find the corresponding fft frequency, and then convert to Hertz:

import wave
import struct
import numpy as np

if __name__ == '__main__':
    data_size = 40000
    fname = "test.wav"
    frate = 11025.0
    wav_file = wave.open(fname, 'r')
    data = wav_file.readframes(data_size)
    wav_file.close()
    data = struct.unpack('{n}h'.format(n=data_size), data)
    data = np.array(data)

    w = np.fft.fft(data)
    freqs = np.fft.fftfreq(len(w))
    print(freqs.min(), freqs.max())
    # (-0.5, 0.499975)

    # Find the peak in the coefficients
    idx = np.argmax(np.abs(w))
    freq = freqs[idx]
    freq_in_hertz = abs(freq * frate)
    print(freq_in_hertz)
    # 439.8975
19
  • @~unutbu:But can I get the frequency values in Hertz?I want to make wav files.
    – ria
    Sep 12, 2010 at 18:29
  • 1
    @PavelShvechikov: Oops, yes. You are absolutely right. Thanks for the correction.
    – unutbu
    Nov 28, 2014 at 13:27
  • 1
    I found it. Basically my data is 2 channel data but your code may not working for me. Jun 14, 2016 at 11:49
  • 1
    I made the wav generation script channels to 2 and then with the script I am getting the freq specified in the wav generation script. But when I record the same. I am getting exactly half of the peak frequency value. What may I go wrong. Thanks in advance Jun 14, 2016 at 12:26
  • 1
    @unutbu Sorry for the revival of this topic. Your example gets the frequency of a wav file that is a constant 440Hz. What if my wav file has 20 samples, at different frequencies, how do I extract all the frequencies in the order they appear? I am able to plot the DFT and normalization with np.fft.fft(signal), but that gets me the frequencies and the number of times they were observed, not the actual order.
    – Jax Teller
    Oct 23, 2020 at 9:16
49

Here we deal with the Numpy implementation of the fft.

Frequencies associated with DFT values (in python)

By fft, Fast Fourier Transform, we understand a member of a large family of algorithms that enable the fast computation of the DFT, Discrete Fourier Transform, of an equisampled signal.

A DFT converts an ordered sequence of N complex numbers to an ordered sequence of N complex numbers, with the understanding that both sequences are periodic with period N.

In many cases, you think of

  • a signal x, defined in the time domain, of length N, sampled at a constant interval dt
  • its DFT X (here specifically X = np.fft.fft(x)), whose elements are sampled on the frequency axis with a sample rate .

Some definition

  • the period (aka duration²) of the signal x, sampled at dt with N samples is is

     T = dt*N
    
  • the fundamental frequencies (in Hz and in rad/s) of X, your DFT are

     df = 1/T
     dω = 2*pi/T # =df*2*pi
    
  • the top frequency is the Nyquist frequency

     ny = dω*N/2
    

    (NB: the Nyquist frequency is not dω*N

The frequencies associated with a particular element in the DFT

The frequencies corresponding to the elements in X = np.fft.fft(x) for a given index 0<=n<N can be computed as follows:

def rad_on_s(n, N, dω):
    return dω*n if n<N/2 else dω*(n-N)

or in a single sweep

ω = np.array([dω*n if n<N/2 else dω*(n-N) for n in range(N)])

if you prefer to consider frequencies in Hz, s/ω/f/

f = np.array([df*n if n<N/2 else df*(n-N) for n in range(N)])

Using those frequencies

If you want to modify the original signal x -> y applying an operator in the frequency domain in the form of a function of frequency only, the way to go is computing the ω's and

Y = X*f(ω)
y = ifft(Y)

Introducing np.fft.fftfreq

Of course numpy has a convenience function np.fft.fftfreq that returns dimensionless frequencies rather than dimensional ones but it's as easy as

f = np.fft.fftfreq(N)*N*df
ω = np.fft.fftfreq(N)*N*dω

Because df = 1/T and T = N/sps (sps being the number of samples per second) one can also write

f = np.fft.fftfreq(N)*sps

Notes

  1. Dual to the sampling interval dt there is the sampling rate, sr, or how many samples are taken during a unit of time; of course dt=1/sr and sr=1/dt.
  2. Speaking of a duration, even if it is rather common, hides the fundamental idea of periodicity.
  3. The concept of Nyquist frequency is clearly exposed in any textbook dealing with the analysis of time signals, and also in the linked Wikipedia article. Does it suffice to say that information cannot be created?
5
  • why there need the condition: n<N/2?
    – Doerthous
    Aug 29, 2021 at 5:31
  • why ny = dw*N/2
    – Doerthous
    Aug 29, 2021 at 7:39
  • @Doerthous Because X is periodic in the frequency domain with period N, you cannot tell if a component of X must be assigned to a time domain signal with frequency n*dw or a frequency (n+k*N)*dw. Nyquist, predating the FFT algorithm, in the context of information theory noted this behavior and introduced the idea of a limit frequency for which the SAMPLED reconstructed signal is identical to a signal reconstructed taking into account higher frequencies (the keyword here is SAMPLED). A decent numerical analysis textbook (or even Wikipedia) will provide you relevant context.
    – gboffi
    Aug 29, 2021 at 8:45
  • After many googled, I'm beginning to catch on. Let fs be the sample frequency which is 1/dt, components whose frequency f >= fs/2 are aliasing to negative frequencies, that is f = df*n < fs/2 = 1/(2dt) = N/2T = df*N/2 => n < N/2, is that right?
    – Doerthous
    Aug 29, 2021 at 14:12
  • @gboffi, how to find the harmonic magnitudes and phase of a 2D DFT image in k-space?
    – pluto
    Feb 15, 2022 at 19:42
7

The frequency is just the index of the array. At index n, the frequency is 2πn / the array's length (radians per unit). Consider:

>>> numpy.fft.fft([1,2,1,0,1,2,1,0])
array([ 8.+0.j,  0.+0.j,  0.-4.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+4.j,
        0.+0.j])

the result has nonzero values at indices 0, 2 and 6. There are 8 elements. This means

       2πit/8 × 0       2πit/8 × 2       2πit/8 × 6
    8 e           - 4i e           + 4i e
y ~ ———————————————————————————————————————————————
                          8
5
  • I'm sorry. But I couldn't get it clearly.Can you tell me what are 't' and 'e' above? why did you introduce 'i*t' in the equation 2πn/8, Is there a function in SciPy doing this calculation?
    – ria
    Sep 12, 2010 at 15:48
  • 1
    @ria: e is 2.71828.... See en.wikipedia.org/wiki/Euler%27s_formula. t is the index of the original array, e.g. t=0 -> 1, t=1 -> 2, t=2 -> 1, etc. Basically, if you want to get the frequency, they are just 0/8, 1/8, 2/8, ..., 7/8.
    – kennytm
    Sep 12, 2010 at 16:05
  • @ KennyTM:I see,the exponent 'e'.I understood.
    – ria
    Sep 12, 2010 at 18:30
  • 2
    This is incorrect – the output of the FFT is not in normal frequency order. See docs.scipy.org/doc/numpy-1.10.0/reference/…
    – jakevdp
    Nov 22, 2015 at 11:51
  • Non-version specific link to numpy FFT documentation: numpy.org/doc/stable/reference/routines.fft.html
    – Chris
    Apr 27, 2020 at 4:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.