Оптимизация числового расчета для аудио программы DSP

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 возвращает массив кортежей, где первое значение - это частота, а второе - комплексное число, представляющее то, что я не совсем понимаю, но в соответствии со страницей, на которой я следовал, чтобы узнать это, это не важно. Вот где я научился этому: 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.])