I used fft function in numpy which resulted in a complex array. How to get the exact frequency values?
- 14,536
- 19
- 67
- 128
- 525
- 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
- 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 -
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
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 atdtwithNsamples is isT = dt*Nthe fundamental frequencies (in Hz and in rad/s) of
X, your DFT aredf = 1/T dw = 2*pi/T # =df*2*pithe 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
- 19,456
- 5
- 52
- 81
-
-
-
@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
- 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
-
-
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