51

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

Ashish Gupta
  • 14,536
  • 19
  • 67
  • 128
ria
  • 525
  • 1
  • 5
  • 6

3 Answers3

73

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
unutbu
  • 777,569
  • 165
  • 1,697
  • 1,613
  • @~unutbu:But can I get the frequency values in Hertz?I want to make wav files. – ria Sep 12 '10 at 18:29
  • @ria: `freq` as defined by `np.fft.fft` is unitless, and always in the interval [-1/2, 1/2]. I believe to convert to Hertz, you multiply by the frame rate and take the absolute value. – unutbu Sep 12 '10 at 19:32
  • No sir, my aim is to read a wav file , perform fft on it,and with each frequency value i have to create a new wav file.i have done till performing fft.this is the output i got: wavefile: Reading fmt chunk Reading data chunk (8000, array([128, 128, 128, ..., 128, 128, 128], dtype=uint8)) numpy array: [128 128 128 ..., 128 128 128] – ria Sep 12 '10 at 20:07
  • fft on numpy: [ 12535945.00000000 +0.j -30797.74496367 +6531.22295858j -26330.14948055-11865.08322966j ..., 34265.08792783+31937.15794965j -26330.14948055+11865.08322966j -30797.74496367 -6531.22295858j] now to write it on a wav file i guess i need the frequency values in hertz.or is there any other way? – ria Sep 12 '10 at 20:08
  • hi,please give me a reply.Because I need to do this ASAP.Is it possible to do this task at least? – ria Sep 13 '10 at 04:51
  • Thanks a lot! That's really helps me. – Mayli Jun 08 '12 at 06:20
  • `idx=np.argmax(np.abs(w)**2)` -- why to square an array of absolute values for finding only index ? you could get same result without squaring. – xolodec Nov 28 '14 at 08:47
  • `np.abs(w)` is an array which contains only real numbers as an absolute value of a complex. So there is no point in squaring them to obtain the maximum of this array. – xolodec Nov 28 '14 at 13:18
  • 1
    @PavelShvechikov: Oops, yes. You are absolutely right. Thanks for the correction. – unutbu Nov 28 '14 at 13:27
  • @unutbu : I will record the tones using arecord or sox's rec. And I saw ur code I was little confused. what is data_size and frate. Can I use ur code for my thing. I also want to find the peak frequency. When I ran the script I got `Traceback (most recent call last): File "aqu.py", line 12, in data = struct.unpack('{n}h'.format(n=data_size), data) struct.error: unpack requires a string argument of length 80000` – optimus prime Jun 14 '16 at 10:53
  • 1
    I found it. Basically my data is 2 channel data but your code may not working for me. – optimus prime Jun 14 '16 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 – optimus prime Jun 14 '16 at 12:26
  • It is working fine for my single channel recorded files. But not for 2 channel recorded files. What I may missed here – optimus prime Jun 14 '16 at 12:37
  • How can you get the frequency values for every fftfreq values instead of the single signal fftfreq value? – Ugur Aug 08 '16 at 21:19
  • Shouldn't the wav file size increase if I switch to `nchannels=2`? – Furqan Hashim Nov 16 '18 at 18:16
  • @FurqanHashim: In the code above, if you change `nchannels=1` to `nchannels=2` then the resultant file will have 2 channels but half as many frames. So the file size does not increase. If `x` and `y` are two wave files (each returned by `wave.open`), you can compare their channels and frames by inspecting `[f.getparams() for f in [x, y]]`. You can use this to confirm the relationship between channels and frames. – unutbu Nov 16 '18 at 18:34
  • @unutbu thanks for the clarification. Just to be sure I assume that `frate` is the number of observations per second? For e.g. in case of stock price, if I have data having a frequency of a second (1 observation per second) `frate` would be 1? And `sampwidth=2` is sampling frequency of 44100 Hz? – Furqan Hashim Nov 16 '18 at 19:01
  • @FurqanHashim: The `frate` is the number of frames per second. So if one stock price corresponds to one frame, then yes, the `frate` would be 1. [The `sampwidth`](https://docs.python.org/3/library/wave.html#wave.Wave_read.getsampwidth) is the number of bytes used to encode one piece of data (or, in your terminology, perhaps one "observation"). Note that there are `nchannel` observations per frame. – unutbu Nov 16 '18 at 20:13
  • How can implement a bandpass filter to this and then obtain the peak frequency? – iamvinitk May 05 '19 at 05:53
  • 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 '20 at 09:16
46

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 a list of N complex numbers to a list of N complex numbers, with the understanding that both lists are periodic with period N.

Here we deal with the numpy implementation of the fft.

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 dw.

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
    dw = 2*pi/T # =df*2*pi
    
  • the top frequency is the Nyquist frequency

    ny = dw*N/2
    

    (and it's not dw*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, dw):
    return dw*n if n<N/2 else dw*(n-N)

or in a single sweep

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

if you prefer to consider frequencies in Hz, s/w/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 w's and

Y = X*f(w)
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
w = np.fft.fftfreq(N)*N*dw

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
gboffi
  • 19,456
  • 5
  • 52
  • 81
  • why there need the condition: n – Doerthous Aug 29 '21 at 05:31
  • why ny = dw*N/2 – Doerthous Aug 29 '21 at 07: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 '21 at 08: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 '21 at 14:12
  • @gboffi, how to find the harmonic magnitudes and phase of a 2D DFT image in k-space? – pluto Feb 15 '22 at 19:42
6

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
kennytm
  • 491,404
  • 99
  • 1,053
  • 989
  • 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 '10 at 15:48
  • 1
    @ria: e is 2.71828.... See http://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 '10 at 16:05
  • @ KennyTM:I see,the exponent 'e'.I understood. – ria Sep 12 '10 at 18:30
  • 2
    This is incorrect – the output of the FFT is not in normal frequency order. See http://docs.scipy.org/doc/numpy-1.10.0/reference/routines.fft.html#implementation-details – jakevdp Nov 22 '15 at 11:51
  • Non-version specific link to numpy FFT documentation: https://numpy.org/doc/stable/reference/routines.fft.html – Chris Apr 27 '20 at 04:35