优化音频 DSP 程序的 numpy 计算

Posted 2023-02-25

【中文标题】优化音频 DSP 程序的 numpy 计算

【英文标题】：Optimizing a numpy calculation for an audio DSP program

【发布时间】：2018-05-08 04:22:21

【问题描述】：

我是一名音乐家，我正在编写一个读取 .wav 文件的 Python 脚本，使用快速傅立叶变换将其转换为一堆正弦波，然后将这些正弦波调谐到最接近的谐波频率。如果所有这些听起来都是胡言乱语，那没关系，我的问题无需任何音乐知识即可回答。

当我在一个相当长的 .wav 文件上运行我的脚本时，需要几个小时来处理脚本的以下部分：

filtered_data_fft = np.zeros(data_fft.size)
for f in data_fft:
    if f > 1:
        valid_frequency = (np.abs(valid_frequencies - i)).argmin()
        filtered_data_fft[valid_frequency] += data_fft[i]
    i += 1

以 fft 结尾的两个数组都是索引对应于频率的数组，valid_frequencies 数组是对应于所述索引的频率列表。我最初并没有对所有东西都使用 numpy 数组，而且运行时间太长，以至于我无法在合理的时间内处理一个简短的声音文件，但是使用 numpy 会快很多。谁能想出比这更快的方法？我会把完整的脚本放在下面。

此外，关于将复数转换为实数的两个已知警告会丢弃复数，但我认为它们不是问题。 FFT 返回一个元组数组，其中第一个值是频率，第二个值是一个复数，表示我不太了解的东西，但根据我学习这一点的页面，这并不重要。这是我学到这些东西的地方：https://pythonforengineers.com/audio-and-digital-signal-processingdsp-in-python/

诚然，我并不完全理解我在这里所做的很多 DSP 工作，所以如果我在某些方面有严重错误，请告诉我！我只是想用一种有趣的方式将噪音变成音乐，用于我正在进行的项目。

这是我正在测试的音频示例： https://my.mixtape.moe/iltlos.wav （重命名为missile.wav）

这是完整的脚本（更新为正确）：

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


# import data from wave
wav_file = wave.open("missile.wav", 'r')
num_samples = wav_file.getnframes()
sampling_rate = wav_file.getframerate() / 2
data = wav_file.readframes(num_samples)
wav_file.close()

data = struct.unpack('nh'.format(n=num_samples), data)
data = np.array(data)

# fast fourier transform makes an array of the frequencies of sine waves that comprise the sound
data_fft = np.fft.rfft(data)


# generate list of ratios that can be used for tuning (not octave reduced)
MAX_HARMONIC = 5
valid_ratios = []
for i in range(1, MAX_HARMONIC + 1):
    for j in range(1, MAX_HARMONIC + 1):
        if i % 2 != 0 and j % 2 != 0:
            valid_ratios.append(i/float(j))
            valid_ratios.append(j/float(i))


# remove dupes
valid_ratios = list(set(valid_ratios))


# find all the frequencies with the valid ratios
valid_frequencies = []
multiple = 2
while(multiple < num_samples / 2):
    multiple *= 2

    for ratio in valid_ratios:
        frequency = ratio * multiple

        if frequency < num_samples / 2:
            valid_frequencies.append(frequency)



# remove dupes and sort and turn into a numpy array
valid_frequencies = np.sort(np.array(list(set(valid_frequencies))))


# bin the data_fft into the nearest valid frequency
valid_frequencies = valid_frequencies.astype(int)
boundaries = np.concatenate([[0], np.round(np.sqrt(0.25 + valid_frequencies[:-1] * valid_frequencies[1:])).astype(int)])
select = np.abs(data_fft) > 1
filtered_data_fft = np.zeros_like(data_fft)
filtered_data_fft[valid_frequencies] = np.add.reduceat(np.where(select, data_fft, 0), boundaries)


# do the inverse fourier transform to get a sound wave back
recovered_signal = np.fft.irfft(filtered_data_fft)

# write sound wave to wave file
comptype="NONE"
compname="not compressed"
nchannels=1
sampwidth=2

wav_file=wave.open("missile_output.wav", 'w')
wav_file.setparams((nchannels, sampwidth, int(sampling_rate), num_samples, comptype, compname))

for s in recovered_signal:
    wav_file.writeframes(struct.pack('h', s))

wav_file.close()

【参考方案1】：

关于脚本的几点说明：

(1) 由于您使用的是rfft，匹配的逆将是irfft 而不是ifft

(2) 就目前而言，脚本接受除了0 之外的每个 频率为有效（因为1 包含在valid_ratios 中

(3) 给定频率的复数包含该“正弦波”的幅度和相位（偏移）。我假设您想根据幅度进行过滤。为此，您必须取复数的绝对值，即np.abs(f) > 1 等。

(4) 一旦你有了一组好的有效频率，你就可以进行如下操作。我同意@Mateen Ulhaq 使用几何中点的建议。

boundaries = np.concatenate([[0], np.round(np.sqrt(0.25 + valid_frequencies[:-1] * valid_frequencies[1:])).astype(int)])
select = np.abs(data_fft) > 1
filtered_data_fft = np.zeros_like(data_fft)
filtered_data_fft[valid_frequencies] = np.add.reduceat(np.where(select, data_fft, 0), boundaries)

【讨论】：

我在你写的最后一行的开头收到了一个错误过滤数据_fft[valid_frequencies] 因为valid_frequencies不是一个布尔数组，我如何将频率列表转换为一个布尔列表，其中所有真值在有效频率的索引号处？

这是查找频率的最新方法，我完全错了paste.ee/p/TgDju

@halbe 你能试试valid_frequencies = valid_frequencies.astype(int)吗？那应该可以。

喜欢这个？ paste.ee/p/LIt6P 我现在收到这个错误：IndexError: index 74316 out-of-bounds in add.reduceat [0, 74305)

@halbe 这是因为 rfft 的输出在 Nyquist 处被剪切（这恰好是一个好主意，无论如何），您可以检查它是输入大小的一半。在计算valid_frequencies 的块中，将num_samples 的两个出现替换为num_samples // 2。

【参考方案2】：

您正在尝试对数据进行分类或数字化。首先定义您的决策边界：

valid_frequencies = np.sort(valid_frequencies)
b = valid_frequencies
b = (b[1:] + b[:-1]) / 2
bins = np.concatenate(([0], b, [MAX_FREQ]))

虽然使用几何平均值而不是算术平均值可能会更成功。 （频率分析通常更多是基于日志的。）

b = np.sqrt(b[1:] * b[:-1])

现在您只需将数据数字化，然后计算各种索引的出现次数：

hist = np.bincount(np.digitize(data_fft, bins))[1:]

也许更快的是：

hist = np.histogram(data_fft, bins=bins)[0]

最后，我们将这些嵌入到正确的索引中：

filtered_data_fft = np.zeros_like(data_fft)
filtered_data_fft[valid_frequencies] = hist

---

编辑：例如，

>>> data_fft = np.array([1.1, 2.2, 3.3, 4.4, 5.5, 6.6, 7.7, 8.8, 9.9])
>>> valid_frequencies = np.sort([2, 4])

>>> b = valid_frequencies
>>> b = (b[1:] + b[:-1]) / 2
>>> bins = np.concatenate(([0.0], b, [10.0]))
array([ 0.,  3., 10.])

>>> hist = np.bincount(np.digitize(data_fft, bins))[1:]
array([2, 7])

>>> hist = np.histogram(data_fft, bins=bins)[0]
array([2, 7])

>>> filtered_data_fft = np.zeros(11)
>>> filtered_data_fft[valid_frequencies] = hist
array([0., 0., 2., 0., 7., 0., 0., 0., 0., 0., 0.])

【讨论】：

它处理得更快，但输出全为零，我错过了什么吗？ paste.ee/p/TbNoT

@halbe 我想你需要小心np.bincount() 在开头输出一个额外的0。我在上面用[1:] 把它切掉了。还要确保num_samples 确实是您的最大频率。或者考虑使用np.histogram，因为这似乎可以更好地处理这些边缘情况。

filtered_data_fft 需要和data_fft 一样大小，它的值是volumes，频率是index 位置。所以 440hz 在索引 440 处，其值是 440hz 正弦波的音量，但在 valid_frequencies 数组中的任意位置都有一个值 440。

@halbe 这比我想象的要复杂。无论如何，我已经更新了答案。

在您的示例中，过滤后的数组应该是 [0., 6.6, 0., 42.9, 0., 0., 0., 0., 0., 0.]) 因为前两个或者三个频率（四舍五入）应该将它们的体积添加到第二个 bin（没有频率 0），其余的应该添加到第四个 bin，因为它是最接近的一个。现在它加起来的是被分箱的频率数量，而不是这些频率的音量。否则它非常接近！我怎样才能对这些值进行分类？