如何在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

在这里，我们讨论的是Numpy中fft的实现。

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

所谓fft，即快速傅里叶变换，是一类算法的统称，这些算法可以快速计算一个等间隔信号的离散傅里叶变换（DFT）。

DFT的作用是将一个有序的个复数序列转换成另一个有序的个复数序列，并且这两个序列都是周期性的，周期为。

在很多情况下，你可以把它想象成：

• 一个信号x，在时间域中定义，长度为，以固定的时间间隔dt进行采样，¹
• 它的DFTX（这里具体指X = np.fft.fft(x)），其元素在频率轴上以采样率dω进行采样。

一些定义

• 信号x的周期（也叫持续时间²），以dt采样，样本数为N，可以表示为

   T = dt*N
  

• DFTX的基本频率（以Hz和rad/s为单位）是

   df = 1/T
   dω = 2*pi/T # =df*2*pi
  

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

   ny = dω*N/2
  

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

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

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

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

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)

---

提问者想知道如何找到以赫兹为单位的频率。 我认为公式是 频率 (Hz) = abs(fft_freq * 帧率)。

下面是一些演示这个过程的代码。

首先，我们创建一个440赫兹的波形文件：

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