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

Question

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

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询问如何以赫兹为单位找到频率.我相信这个公式是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 2

与DFT值相关的频率(在python中)

通过fft,快速傅立叶变换,我们了解了一大类算法的成员,这些算法能够快速计算等采样信号的DFT离散傅立叶变换.

一个DFT列表转换ň复数到列表ñ复数,但有一项谅解,这两个名单是周期性的,周期ñ.

在这里我们处理fft的numpy实现.

在很多情况下,你会想到

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

一些定义

• 的信号的周期(也称为持续时间)x,在取样dt与N样品是被

  T = dt*N
  

• X您的DFT 的基频(以Hz和rad/s为单位)

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

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

  ny = dw*N/2
  

  (而且不是dw*N)

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

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

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

或者一次扫描

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

如果您更愿意考虑以Hz为单位的频率, 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

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