优化音频 DSP 程序的 numpy 计算

Question

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

当我运行我的脚本在一个相当长的 .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 returns 一个元组数组，其中第一个值是一个频率，第二个是一个复数，代表我不太了解的东西，但根据我遵循的页面了解这一点并不重要。这是我学到这些东西的地方：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

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

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

Answer 2

关于您的脚本的几点说明：

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

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

(3) 给定频率的复数包含 "sine wave" 的振幅和相位（位移）。我假设您想根据振幅进行过滤。为此，您必须取复数的绝对值，即 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)

优化音频 DSP 程序的 numpy 计算

Optimizing a numpy calculation for an audio DSP program

python

arrays

optimization

signal-processing

numpy