如何在python中提取与fft值相关的频率

Question

我在 numpy 中使用了 fft 函数，这导致了一个复杂的数组。如何得到准确的频率值？

Answer 1

频率只是数组的索引。在索引 n 处，频率为 2πn / 数组的长度（每单位的弧度）。考虑：

>>> 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])

结果在索引 0、2 和 6 处具有非零值。有 8 个元素。这意味着

       2πit/8 × 0       2πit/8 × 2       2πit/8 × 6
    8 e           - 4i e           + 4i e
y ~ ———————————————————————————————————————————————
                          8

Answer 2

np.fft.fftfreq 告诉您与系数相关的频率：

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)

---

OP 询问如何以赫兹为单位找到频率。 我相信公式是frequency (Hz) = abs(fft_freq * frame_rate)。

这里有一些代码可以证明这一点。

首先，我们制作一个 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()

这将创建文件test.wav。 现在我们读入数据，对其进行 FFT，找到具有最大功率的系数， 并找到对应的fft频率，然后转换为赫兹：

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

Answer 3

与 DFT 值相关的频率（在 python 中）

fft，快速傅里叶变换，我们了解一个大型算法家族的成员，这些算法能够快速计算等采样的 DFT（离散傅里叶变换）信号。

A DFT 将 N 个复数列表转换为 N 个复数列表，同时理解这两个列表都是周期性的，周期为 N .

这里我们处理fft的numpy实现。

很多时候你会想到

• 在长度为 N 的时域中定义的信号 x，在 a 恒定间隔dt,
• 它的 DFT X（这里特指X = np.fft.fft(x)），其元素 以 dw 的采样率在频率轴上进行采样。

一些定义

• 信号x的周期（又名持续时间），在dt采样，N采样是

  T = dt*N
  

• X 的基本频率（以 Hz 和 rad/s 为单位），您的 DFT 是

  df = 1/T
  dw = 2*pi/T # =df*2*pi
  

• 最高频率是Nyquist frequency

  ny = dw*N/2
  

  （而且不是dw*N）

与 DFT 中特定元素相关的频率

对于给定索引0<=n<N，与X = np.fft.fft(x) 中的元素对应的频率可以计算如下：

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

或单次扫描

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

如果您更喜欢以赫兹为单位考虑频率，s/w/f/

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

使用这些频率

如果您想修改原始信号x -> y 在频域中仅以频率函数的形式应用算子，则可以计算w 和

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

介绍np.fft.fftfreq

当然numpy 有一个方便的函数np.fft.fftfreq，它返回无量纲频率而不是维度频率，但它就像

一样简单

f = np.fft.fftfreq(N)*N*df
w = np.fft.fftfreq(N)*N*dw

因为df = 1/T 和T = N/sps（sps 是每秒的样本数）也可以这样写

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