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

- 14,869
- 20
- 75
- 134

- 555
- 1
- 5
- 6
3 Answers
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

- 842,883
- 184
- 1,785
- 1,677
-
@~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
-
`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 -
1I 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
-
1I 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
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 dω.
Some definition
the period (aka duration²) of the signal
x
, sampled atdt
withN
samples is isT = dt*N
the fundamental frequencies (in Hz and in rad/s) of
X
, your DFT aredf = 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
- 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.
- Speaking of a duration, even if it is rather common, hides the fundamental idea of periodicity.
- 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?

- 22,939
- 8
- 54
- 85
-
-
-
@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
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

- 510,854
- 105
- 1,084
- 1,005
-
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
-
-
2This 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