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

Question

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

我使用Audacity绘制并验证来自相同440Hz wav文件的光谱，我得到了这个，这表明440Hz确实是主导频率：（https://i.imgur.com/2UImEkR.png）

要使用python执行此操作，我使用PyAudio库并引用this blog。到目前为止，我使用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.imgur.com/zsAXME5.png）

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

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

Answer 1

Forewords

正如xdurch0指出的那样，你正在阅读一种索引而不是频率。如果您要自己进行所有计算，则需要在绘图之前计算自己的频率向量，如果您想获得一致的结果。阅读本answer可以帮助您解决问题。

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

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

或（全机）：

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

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

MCVE

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

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

基本上：

• Read the wav file并获得采样率（此处为44.1kHz）;
• 计算Power Spectrum Density和频率;
• 然后用matplotlib显示它。

这输出：

Find Peak

然后我们可以使用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]]

Handling multiple channels

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

Wav file可以容纳多个频道，主要有单声道或立体声文件（最大2**16 - 1频道）。由于多个通道文件（stereo sample）而出现下划线的问题。

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

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

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

它也适用于转置data.T但我们需要回调转置结果。

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

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

Real application

这种方法适用于真正的分块样本，只要块有足够的数据（参见Nyquist-Shannon sampling theorem）。然后，数据是信号的子样本（块），并且速率保持不变，因为它在过程中不会改变。

有大小的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进行后处理，以便更容易识别。