从 PyAudio 接收的数据的 FFT 给出错误的频率

Question

我的主要任务是实时识别麦克风中的人类嗡嗡声。作为识别一般信号的第一步，我对手机上的应用程序生成的 440 Hz 信号进行了 5 秒的记录，并尝试检测相同的频率。

我使用 Audacity 绘制并验证了同一个 440Hz wav 文件的频谱，我得到了这个，这表明 440Hz 确实是主频率：( https://i.stack.imgur.com/c3DWD.png )

为了使用 python 执行此操作，我使用PyAudio库并参考此博客。到目前为止，我使用 wav 文件运行的代码是这样的：

"""PyAudio Example: Play a WAVE file."""

import pyaudio
import wave
import sys
import struct
import numpy as np
import matplotlib.pyplot as plt

CHUNK = 1024

if len(sys.argv) < 2:
    print("Plays a wave file.\n\nUsage: %s filename.wav" % sys.argv[0])
    sys.exit(-1)

wf = wave.open(sys.argv[1], 'rb')

p = pyaudio.PyAudio()
stream = p.open(format=p.get_format_from_width(wf.getsampwidth()),
                channels=wf.getnchannels(),
                rate=wf.getframerate(),
                output=True)

data = wf.readframes(CHUNK)

i = 0
while data != '':
    i += 1
    data_unpacked = struct.unpack('{n}h'.format(n= len(data)/2 ), data) 
    data_np = np.array(data_unpacked) 
    data_fft = np.fft.fft(data_np)
    data_freq = np.abs(data_fft)/len(data_fft) # Dividing by length to normalize the amplitude as per https://www.mathworks.com/matlabcentral/answers/162846-amplitude-of-signal-after-fft-operation
    print("Chunk: {} max_freq: {}".format(i,np.argmax(data_freq)))

    fig = plt.figure()
    ax = fig.add_subplot(1,1,1)
    ax.plot(data_freq)
    ax.set_xscale('log')
    plt.show()

    stream.write(data)
    data = wf.readframes(CHUNK)

stream.stop_stream()
stream.close()

p.terminate()

在输出中，我得到所有块的最大频率为 10，其中一个图的示例是：( https://i.stack.imgur.com/2e3wR.png )

我原本预计所有块的这个值都是 440，而不是 10。我承认我对 FFT 的理论知之甚少，并且感谢任何帮助我解决这个问题的帮助。

编辑：采样率为 44100。通道数为 2，样本宽度也为 2。

Answer 1

前言

正如所xdurch0指出的，您正在阅读一种索引而不是频率。如果您要自己进行所有计算，如果您想获得一致的结果，则需要在绘图之前计算自己的频率向量。阅读这个答案可能会帮助您找到解决方案。

FFT（半平面）的频率向量为：

 f = np.linspace(0, rate/2, N_fft/2)

或者（完整平面）：

 f = np.linspace(-rate/2, rate/2, N_fft)

另一方面，我们可以将大部分工作委托给优秀的scipy.signal工具箱，该工具箱旨在解决此类问题（以及更多问题）。

MCVE

使用scipy包可以直接获得WAV具有单一频率的简单文件的所需结果（源）：

import numpy as np
from scipy import signal
from scipy.io import wavfile
import matplotlib.pyplot as plt

# Read the file (rate and data):
rate, data = wavfile.read('tone.wav') # See source

# Compute PSD:
f, P = signal.periodogram(data, rate) # Frequencies and PSD

# Display PSD:
fig, axe = plt.subplots()
axe.semilogy(f, P)
axe.set_xlim([0,500])
axe.set_ylim([1e-8, 1e10])
axe.set_xlabel(r'Frequency, $\nu$ $[\mathrm{Hz}]$')
axe.set_ylabel(r'PSD, $P$ $[\mathrm{AU^2Hz}^{-1}]$')
axe.set_title('Periodogram')
axe.grid(which='both')

基本上：

• 读取wav文件并获取采样率（此处44.1kHz）；
• 计算功率谱密度和频率；
• 然后用 来显示它matplotlib。

这输出：

寻找峰值

然后我们可以使用 找到第一个最高峰值的频率（P>1e-2，该标准需要调整）find_peaks：

idx = signal.find_peaks(P, height=1e-2)[0][0]
f[idx] # 440.0 Hz

把所有这些放在一起，可以归结为：

def freq(filename, setup={'height': 1e-2}):
    rate, data = wavfile.read(filename)
    f, P = signal.periodogram(data, rate)
    return f[signal.find_peaks(P, **setup)[0][0]]

处理多个通道

我用我的 wav 文件尝试了这段代码，并得到了 axe.semilogy(f, Pxx_den) 行的错误，如下所示： ValueError: x 和 y 必须具有相同的第一维度。我检查了形状，f 有 (2,)，而 Pxx_den 有 (220160,2)。此外，Pxx_den 数组似乎只有全零。

Wav 文件可以容纳多个通道，主要有单声道或立体声文件（最大2**16 - 1通道数）。您下划线的问题是由于多通道文件（立体声样本）而发生的。

rate, data = wavfile.read('aaaah.wav') # Shape: (46447, 2), Rate: 48 kHz

它没有很好的记录，但该方法signal.periodogram也在矩阵上执行，并且其输入与输出不直接一致wavfile.read（默认情况下它们在不同的轴上执行）。axis因此，在执行 PSD 时，我们需要仔细定向尺寸（使用开关）：

f, P = signal.periodogram(data, rate, axis=0, detrend='linear')

它也适用于转置data.T，但我们需要对结果进行反向转置。

指定轴解决了问题：频率向量是正确的，并且 PSD 在任何地方都不为空（在对axis=1长度为 的执行之前2，在您的情况下，它对我们想要的 2 个样本信号执行 220160 PSD）。

该detrend开关确保信号具有零均值并消除其线性趋势。

实际应用

这种方法应该适用于真正的分块样本，前提是块包含足够的数据（请参阅奈奎斯特-香农采样定理）。然后数据是信号（块）的子样本，并且速率保持恒定，因为它在过程中不会改变。

拥有大小的块2**10似乎很有效，我们可以从中识别特定的频率：

f, P = signal.periodogram(data[:2**10,:], rate, axis=0, detrend='linear') # Shapes: (513,) (513, 2)
idx0 = signal.find_peaks(P[:,0], threshold=0.01, distance=50)[0] # Peaks: [46.875, 2625., 13312.5, 16921.875] Hz

fig, axe = plt.subplots(2, 1, sharex=True, sharey=True)
axe[0].loglog(f, P[:,0])
axe[0].loglog(f[idx0], P[idx0,0], '.')
# [...]

此时，最棘手的部分是微调find-peaks方法以捕获所需的频率。您可能需要考虑对信号进行预过滤或对 PSD 进行后处理，以便更容易识别。