优化音频DSP程序的numpy计算

我是一名音乐家,我正在编写一个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\u frequencies阵列是对应于所述索引的频率列表。最初,我并没有在所有方面都使用numpy阵列,而且运行时间太长,以至于我无法在合理的时间内处理一个简短的声音文件,但使用numpy速度要快得多。有人能想出一个比这更快的方法吗?我将把完整的脚本放在下面。

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

诚然,我并不完全理解我在这里做的很多DSP的东西,所以如果我在某些方面大错特错,请告诉我!我只是想为我正在进行的一个项目找到一种有趣的方式,将噪音转化为音乐。

以下是我正在测试的音频示例: https://my.mixtape.moe/iltlos.wav (将其重命名为导弹.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()