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

Question

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

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

要使用 python 执行此操作，我使用 PyAudio library and refer 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

前言

正如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 包可以直接获得具有单一频率 (source) 的简单 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')

基本上：

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

这输出：

找到峰值

然后我们可以使用find_peaks:

找到第一个最高峰的频率（P>1e-2，这个标准需要调整）

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]]

处理多个渠道

I tried this code with my wav file, and got the error for the line axe.semilogy(f, Pxx_den) as follows : ValueError: x and y must have same first dimension. I checked the shapes and f has (2,) while Pxx_den has (220160,2). Also, the Pxx_den array seems to have all zeros only.

Wav file can hold multiple channels, mainly there are mono or stereo files (max. 2**16 - 1 channels). The problem you underlined occurs because of multiple channels file (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 不是到处都是空的（在它对长度为 2 的 axis=1 上执行之前，在你的情况下它执行了 220160 PSD在 2 样本信号上我们想要相反的结果。

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

实际应用

如果块包含足够的数据（参见 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 方法以捕获所需的频率。您可能需要考虑预过滤信号或 post 处理 PSD 以便更容易识别。

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

FFT of data received from PyAudio gives wrong frequency

python

signal-processing

audio-processing

pyaudio

前言

MCVE

找到峰值

处理多个渠道

实际应用