从麦克风产生频谱图

Question

下面我的代码将从麦克风获取输入，如果音频块的平均值超过某个阈值，它将生成音频块的频谱图（长 30 毫秒）。以下是正常对话中生成的频谱图的样子：

据我所见，考虑到音频及其环境，这看起来不像我期望的频谱图。我期待更像以下的东西（转置以节省空间）：

我正在录制的麦克风是我的 Macbook 上的默认麦克风，有什么建议可以解决问题吗？

---

record.py：

import pyaudio
import struct
import math
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt


THRESHOLD = 40 # dB
RATE = 44100
INPUT_BLOCK_TIME = 0.03 # 30 ms
INPUT_FRAMES_PER_BLOCK = int(RATE * INPUT_BLOCK_TIME)

def get_rms(block):
    return np.sqrt(np.mean(np.square(block)))

class AudioHandler(object):
    def __init__(self):
        self.pa = pyaudio.PyAudio()
        self.stream = self.open_mic_stream()
        self.threshold = THRESHOLD
        self.plot_counter = 0

    def stop(self):
        self.stream.close()

    def find_input_device(self):
        device_index = None
        for i in range( self.pa.get_device_count() ):
            devinfo = self.pa.get_device_info_by_index(i)
            print('Device %{}: %{}'.format(i, devinfo['name']))

            for keyword in ['mic','input']:
                if keyword in devinfo['name'].lower():
                    print('Found an input: device {} - {}'.format(i, devinfo['name']))
                    device_index = i
                    return device_index

        if device_index == None:
            print('No preferred input found; using default input device.')

        return device_index

    def open_mic_stream( self ):
        device_index = self.find_input_device()

        stream = self.pa.open(  format = pyaudio.paInt16,
                                channels = 1,
                                rate = RATE,
                                input = True,
                                input_device_index = device_index,
                                frames_per_buffer = INPUT_FRAMES_PER_BLOCK)

        return stream

    def processBlock(self, snd_block):
        f, t, Sxx = signal.spectrogram(snd_block, RATE)
        plt.pcolormesh(t, f, Sxx)
        plt.ylabel('Frequency [Hz]')
        plt.xlabel('Time [sec]')
        plt.savefig('data/spec{}.png'.format(self.plot_counter), bbox_inches='tight')
        self.plot_counter += 1

    def listen(self):
        try:
            raw_block = self.stream.read(INPUT_FRAMES_PER_BLOCK, exception_on_overflow = False)
            count = len(raw_block) / 2
            format = '%dh' % (count)
            snd_block = np.array(struct.unpack(format, raw_block))
        except Exception as e:
            print('Error recording: {}'.format(e))
            return

        amplitude = get_rms(snd_block)
        if amplitude > self.threshold:
            self.processBlock(snd_block)
        else:
            pass

if __name__ == '__main__':
    audio = AudioHandler()
    for i in range(0,100):
        audio.listen()

---

基于 cmets 的编辑：

如果我们将速率限制为 16000 Hz 并为颜色图使用对数刻度，则这是在麦克风附近轻敲的输出：

这在我看来仍然有些奇怪，但似乎也是朝着正确方向迈出的一步。

使用 Sox 并与我的程序生成的频谱图进行比较：

Answer 1

首先，观察您的代码最多可以绘制 100 个频谱图（如果多次调用 processBlock），您只能看到最后一个。你可能想解决这个问题。此外，我假设您知道为什么要使用 30 毫秒的录音。就个人而言，我想不出用笔记本电脑麦克风记录 30 毫秒的实际应用可以提供有趣的见解。这取决于您录制的内容以及触发录制的方式，但这个问题与实际问题无关。

否则代码可以完美运行。只需对 processBlock 函数进行一些小改动，应用一些背景知识，您就可以获得信息丰富且美观的频谱图。

让我们来谈谈实际的频谱图。我将以 SoX 输出作为参考。颜色条注释说它是dBFS¹，这是一个对数度量（dB 是Decibel 的缩写）。所以，让我们先将频谱图转换为 dB：

    f, t, Sxx = signal.spectrogram(snd_block, RATE)   
    dBS = 10 * np.log10(Sxx)  # convert to dB
    plt.pcolormesh(t, f, dBS)

这改进了色阶。现在我们看到了之前隐藏的较高频段的噪声。接下来，让我们解决时间分辨率问题。频谱图将信号分成段（默认长度为 256）并计算每个段的频谱。这意味着我们具有出色的频率分辨率，但时间分辨率非常差，因为只有少数这样的段适合信号窗口（大约 1300 个样本长）。时间和频率分辨率之间总是需要权衡取舍。这与uncertainty principle 有关。因此，让我们通过将信号分成更短的部分来用一些频率分辨率换取时间分辨率：

f, t, Sxx = signal.spectrogram(snd_block, RATE, nperseg=64)

太棒了！现在我们在两个轴上都得到了一个相对平衡的分辨率——但是等等！为什么结果如此像素化？！实际上，这就是短短 30 毫秒时间窗口中的所有信息。 1300 个样本在二维中分布的方式只有这么多。但是，我们可以作弊并使用更高的 FFT 分辨率和重叠段。这使得结果更平滑，尽管它没有提供额外的信息：

f, t, Sxx = signal.spectrogram(snd_block, RATE, nperseg=64, nfft=256, noverlap=60)

看漂亮的光谱干涉图案。 （这些模式取决于所使用的窗口函数，但我们不要在这里详细介绍。请参阅频谱图函数的window 参数来玩这些。）结果看起来不错，但实际上不包含更多信息上一张图片。

为了使结果更加 SoX-lixe 观察 SoX 频谱图在时间轴上相当模糊。您可以通过使用原始的低时间分辨率（长段）来获得这种效果，但让它们重叠以获得平滑度：

f, t, Sxx = signal.spectrogram(snd_block, RATE, noverlap=250)

我个人更喜欢第三种解决方案，但您需要找到自己喜欢的时间/频率折衷方案。

最后，让我们使用更像 SoX 的颜色图：

plt.pcolormesh(t, f, dBS, cmap='inferno')

对以下行的简短评论：

THRESHOLD = 40 # dB

阈值与输入信号的 RMS 进行比较，该输入信号 不是以 dB 为单位，而是以原始幅度单位衡量。

---

¹ 显然 FS 是 full scale 的缩写。 dBFS 表示 dB 测量值是相对于最大范围的。 0 dB 是当前表示中可能的最大信号，因此实际值必须

Answer 2

更新为了使我的答案更清晰，并希望补充 @kazemakase 的出色解释，我发现了三件事希望对您有所帮助：

1. 使用 LogNorm：

   plt.pcolormesh(t, f, Sxx, cmap='RdBu', norm=LogNorm(vmin=Sxx.min(), vmax=Sxx.max()))
   

2. 使用numpy的fromstring方法

结果表明 RMS 计算不适用于此方法，因为数据是长度受限的数据类型，并且溢出变为负数：即 507*507=-5095。

1. 使用 colorbar()，因为当你可以看到比例时，一切都会变得更容易

   plt.colorbar()
   

原答案：

在您的代码中播放 10kHz 频率时，我得到了一个不错的结果，只需进行一些更改：

• 导入 LogNorm

  from matplotlib.colors import LogNorm
  

• 在网格中使用 LogNorm

  plt.pcolormesh(t, f, Sxx, cmap='RdBu', norm=LogNorm(vmin=Sxx.min(), vmax=Sxx.max()))
  

这给了我：

您可能还需要在 savefig 之后调用 plt.close()，并且我认为流读取需要一些工作，因为后来的图像会丢失第一季度的声音。

我还推荐plt.colorbar()，这样你就可以看到它最终使用的规模

更新：看到有人花时间投反对票

这是我的频谱图工作版本的代码。 它捕获五秒钟的音频并将它们写入规范文件和音频文件，以便您进行比较。还有很多需要改进的地方，而且几乎没有优化：我确定它会因为编写音频和规范文件的时间而丢弃块。更好的方法是使用非阻塞回调，我以后可能会这样做

与原始代码的主要区别在于为 numpy 获取正确格式的数据：

np.fromstring(raw_block,dtype=np.int16)

而不是

struct.unpack(format, raw_block)

当我尝试将音频写入文件时，这显然是一个主要问题：

scipy.io.wavfile.write('data/audio{}.wav'.format(self.plot_counter),RATE,snd_block)

这是一段好听的音乐，鼓声很明显：

代码：

import pyaudio
import struct
import math
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
import time
from scipy.io.wavfile import write

THRESHOLD = 0 # dB
RATE = 44100
INPUT_BLOCK_TIME = 1 # 30 ms
INPUT_FRAMES_PER_BLOCK = int(RATE * INPUT_BLOCK_TIME)
INPUT_FRAMES_PER_BLOCK_BUFFER = int(RATE * INPUT_BLOCK_TIME)

def get_rms(block):
    return np.sqrt(np.mean(np.square(block)))

class AudioHandler(object):
    def __init__(self):
        self.pa = pyaudio.PyAudio()
        self.stream = self.open_mic_stream()
        self.threshold = THRESHOLD
        self.plot_counter = 0

    def stop(self):
        self.stream.close()

    def find_input_device(self):
        device_index = None
        for i in range( self.pa.get_device_count() ):
            devinfo = self.pa.get_device_info_by_index(i)
            print('Device %{}: %{}'.format(i, devinfo['name']))

            for keyword in ['mic','input']:
                if keyword in devinfo['name'].lower():
                    print('Found an input: device {} - {}'.format(i, devinfo['name']))
                    device_index = i
                    return device_index

        if device_index == None:
            print('No preferred input found; using default input device.')

        return device_index

    def open_mic_stream( self ):
        device_index = self.find_input_device()

        stream = self.pa.open(  format = self.pa.get_format_from_width(2,False),
                                channels = 1,
                                rate = RATE,
                                input = True,
                                input_device_index = device_index)

        stream.start_stream()
        return stream

    def processBlock(self, snd_block):
        f, t, Sxx = signal.spectrogram(snd_block, RATE)
        zmin = Sxx.min()
        zmax = Sxx.max()
        plt.pcolormesh(t, f, Sxx, cmap='RdBu', norm=LogNorm(vmin=zmin, vmax=zmax))
        plt.ylabel('Frequency [Hz]')
        plt.xlabel('Time [sec]')
        plt.axis([t.min(), t.max(), f.min(), f.max()])
        plt.colorbar()
        plt.savefig('data/spec{}.png'.format(self.plot_counter), bbox_inches='tight')
        plt.close()
        write('data/audio{}.wav'.format(self.plot_counter),RATE,snd_block)
        self.plot_counter += 1

    def listen(self):
        try:
            print "start", self.stream.is_active(), self.stream.is_stopped()
            #raw_block = self.stream.read(INPUT_FRAMES_PER_BLOCK, exception_on_overflow = False)

            total = 0
            t_snd_block = []
            while total < INPUT_FRAMES_PER_BLOCK:
                while self.stream.get_read_available() <= 0:
                  print 'waiting'
                  time.sleep(0.01)
                while self.stream.get_read_available() > 0 and total < INPUT_FRAMES_PER_BLOCK:
                    raw_block = self.stream.read(self.stream.get_read_available(), exception_on_overflow = False)
                    count = len(raw_block) / 2
                    total = total + count
                    print "done", total,count
                    format = '%dh' % (count)
                    t_snd_block.append(np.fromstring(raw_block,dtype=np.int16))
            snd_block = np.hstack(t_snd_block)
        except Exception as e:
            print('Error recording: {}'.format(e))
            return

        self.processBlock(snd_block)

if __name__ == '__main__':
    audio = AudioHandler()
    for i in range(0,5):
        audio.listen()

Answer 3

我认为问题在于您正在尝试制作 30 毫秒音频块的频谱图，该音频块太短以至于您可以将信号视为静止的。
频谱图实际上是 STFT，您也可以在 Scipy documentation 中找到它：

scipy.signal.spectrogram(x, fs=1.0, window=('tukey', 0.25), nperseg=None, noverlap=None, nfft=None, detrend='constant', return_oneside=True, scaling='density', axis=-1, mode='psd')

计算具有连续傅里叶变换的频谱图。

频谱图可用作可视化非平稳信号频率内容随时间变化的一种方式。

在第一个图中，您有四个切片，它们是信号块上四个连续 fft 的结果，带有一些窗口和重叠。第二个图有一个独特的切片，但它取决于您使用的频谱图参数。
关键是你想用那个信号做什么。算法的目的是什么？

Answer 4

我不确定直接在 Python 中工作是否是声音处理的最佳方式，最准确地说是使用 FFT...[在我看来，使用 cython 似乎是使用 python 进行声音处理的义务]

您是否评估了绑定任何外部 FFT 方法的可能性（例如使用 fftw）并继续使用 python 仅将作业分配给外部方法并更新图片结果？

您可以在 python here 中找到一些与优化 FFT 相关的信息，也可以查看 scipy FFT 实现。