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

Question

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

Answer 1

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 询问如何找到以赫兹为单位的频率。
我相信公式是频率（Hz）=abs（fft_freq *帧速率）。

下面是一些代码来演示这一点。

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

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

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

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

DFT 转换 N 的有序序列 em> 复数转换为 N 个复数的有序序列，要理解的是这两个序列都是周期性且周期为 N。

在许多情况下，您会想到

• 在时域中定义的信号x，长度为N，以恒定间隔dt采样，
• 其 DFT X（此处具体为 X = np.fft.fft(x)），其元素在频率轴上以采样率 dω 进行采样嗯>。

一些定义

• 信号x的周期（又名持续时间²），在dt处采样，具有N个样本，是

  <前><代码> T = dt*N

• 基频（以Hz为单位，以rad/s为单位）的X，你的DFT是

  <前><代码> df = 1/T
  dω = 2*pi/T # =df*2*pi

• 最高频率是 奈奎斯特频率

  <前><代码> ny = dω*N/2

  （注意：奈奎斯特频率不是 dω*N）³

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

对于给定索引 0 中的元素相对应的频率=n 可以按如下方式计算：

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

或者在单次扫描中

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

，如果您更愿意考虑以 Hz 为单位的频率，s/ω/f/

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

使用这些频率

如果您想修改原始信号x -> y 仅以频率函数的形式在频域中应用算子，方法是计算 ω 并

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

引入 np.fft .fftfreq

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

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

因为 df = 1/T 和 T = N/sps （sps 是每个样本的样本数） 一样简单第二）还可以写

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

---

注释

1. 与采样间隔dt对偶，有采样率，sr，或者在一个过程中采集了多少个样本时间单位；当然是dt=1/sr和sr=1/dt。
2. 说到持续时间，即使它相当常见，也隐藏了
   周期性的基本思想。
3. 奈奎斯特频率的概念在任何涉及时间信号分析的教科书中以及链接的维基百科文章中都清楚地暴露出来。难道说信息不能被创造就足够了吗？

Answer 3

频率只是数组的索引。在索引 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