优化音频 DSP 程序的 numpy 计算

Question

我是一名音乐家，我正在编写一个读取 .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('{n}h'.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()

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

Answer 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.])