Звучит больно, не правда ли? Что еще говорить о том, что одна и та же нота звучит по-разному на разных музыкальных инструментах. Почему же так? Все дело тут в наличии дополнительных гармоник, создающих уникальный тембр каждого инструмента.
Но нас интересует другой вопрос: как этот уникальный тембр смоделировать на компьютере?
Стандартный алгоритм Карплуса-Стронга

Иллюстрация взята с этого сайта.
Суть алгоритма в следующем:
1) Создаем массив размера N из случайных чисел (N напрямую связана с основной частотой звука).
2) Добавляем к концу этого массива значение, посчитанное по следующей формуле:
$$display$$y(n)=\frac{y(n-N)+y(n-N-1)}{2},$$display$$
где $inline$y$inline$ наш массив.
3) Выполняем пункт 2 необходимое количество раз.
Приступим к написанию кода:
1) Импортируем необходимые библиотеки.
import numpy as npimport scipy.io.wavfile as wave
2) Инициализируем переменные.
frequency = 82.41 # Основная частота сигнала в Гцduration = 1 # Время сигнала в секундахsample_rate = 44100 # Частота дискретизации
3) Создаем шум.
# Частота сигнала, равная frequency, означает, что сигнал должен колеблется за одну секунду frequency раз.# Сигнал за одну секунду колеблется sample_rate/length раз.# Тогда length = sample_rate/frequency.noise = np.random.uniform(-1, 1, int(sample_rate/frequency))
4) Создаем массив для хранения значений и добавляем шум в начале.
samples = np.zeros(int(sample_rate*duration))for i in range(len(noise)): samples[i] = noise[i]
5) Используем формулу.
for i in range(len(noise), len(samples)): # В начале i меньше длины шума, поэтому мы берем значения из шума. # Но потом, когда i больше длины шума, мы уже берем посчитанные нами новые значения. samples[i] = (samples[i-len(noise)]+samples[i-len(noise)-1])/2
6) Нормируем и переводим в нужный тип данных.
samples = samples / np.max(np.abs(samples)) samples = np.int16(samples * 32767)
7) Сохраняем в файл.
wave.write("SoundGuitarString.wav", 44100, samples)
8) Оформим все как функцию. Собственно, вот и весь код.
import numpy as npimport scipy.io.wavfile as wave def GuitarString(frequency, duration=1., sample_rate=44100, toType=False): # Частота сигнала, равная frequency, означает, что сигнал должен колеблеться за одну секунду frequency раз. # Сигнал за одну секунду колеблется sample_rate/length раз. # Тогда length = sample_rate/frequency. noise = np.random.uniform(-1, 1, int(sample_rate/frequency)) # Создаем шум samples = np.zeros(int(sample_rate*duration)) for i in range(len(noise)): samples[i] = noise[i] for i in range(len(noise), len(samples)): # В начале i меньше длины шума, поэтому мы берем значения из шума. # Но потом, когда i больше длины шума, мы уже берем посчитанные нами новые значения. samples[i] = (samples[i-len(noise)]+samples[i-len(noise)-1])/2 if toType: samples = samples / np.max(np.abs(samples)) # Нормируем от -1 до 1 return np.int16(samples * 32767) # Переводим в тип данных int16 else: return samples frequency = 82.41sound = GuitarString(frequency, duration=4, toType=True)wave.write("SoundGuitarString.wav", 44100, sound)
9) Запустим и получим:
Для того, чтобы струна звучала лучше, слегка улучшим формулу:
$$display$$y(n)=0.996\frac{y(n-N)+y(n-N-1)}{2}$$display$$
Открытая шестая струна (82.41 Гц) звучит так:
Открытая первая струна (329.63 Гц) звучит так:
Звучит неплохо, не правда ли?
Можно бесконечно подбирать этот коэффициент и найти среднее между красивым звучанием и длительностью, но лучше сразу перейти к Расширенному алгоритму Карплуса-Стронга.
Немного о Z-преобразовании
Пусть $inline$x$inline$ массив входных значений, а $inline$y$inline$ массив выходных значений. Каждый элемент массива y выражается следующей формулой:
$$display$$y(n)=x(n)+x(n-1).$$display$$
Если индекс выходит за пределы массива, то значение равно 0. То есть $inline$x(0-1)=0$inline$. (Посмотрите прошлый код, там это неявно использовалось).
Эту формулу можно записать в соответствующем Z-преобразовании:
$$display$$H(z)=1+z^{-1}.$$display$$
Если формула такая:
$$display$$y(n)=x(n)+x(n-1)-y(n-1).$$display$$
То есть каждый элемент входного массива зависит от прошлого элемента этого же массива (кроме нулевого элемента, конечно). То соответствующее Z-преобразование выглядит так:
$$display$$H(z)=\frac{1+z^{-1}}{1+z^{-1}}.$$display$$
Обратный процесс: из Z-преобразования получить формулу для каждого элемента. Например,
$$display$$H(z)=\frac{1+z^{-1}}{1-z^{-1}}.$$display$$
$$display$$H(z)=\frac{Y(z)}{X(z)} =\frac{1+z^{-1}}{1-z^{-1}}. $$display$$
$$display$$Y(z)*(1-z^{-1} )=X(z)*(1+z^{-1} ). $$display$$
$$display$$Y(z)*1-Y(z)*z^{-1}= X(z)*1+X(z)*z^{-1}.$$display$$
$$display$$y(n)-y(n-1)=x(n)+x(n-1). $$display$$
$$display$$y(n)=x(n)+x(n-1)+y(n-1).$$display$$
Если кто-то не понял, то формула такая: $inline$Y(z)**z^{-k}=*y(n-k)$inline$, где $inline$$inline$ любое действительное число.
Если надо умножить два Z-преобразования друг на друга, то $inline$z^{-a}*z^{-b}=z^{-a-b}.$inline$
Расширенный алгоритм Карплуса-Стронга

Иллюстрация взята с этого сайта.
Вот быстрое описание каждой функции.
Часть I. Функции, преобразующие начальный шум
1) Pick-direction lowpass filter (Фильтр низких частот) $inline$H_p (z)$inline$.
$$display$$H_p (z)=\frac{1-p}{1-pz^{-1}},p [0,1).$$display$$
Соответствующая формула:
$$display$$y(n)=(1-p)x(n)+py(n-1).$$display$$
Код:
buffer = np.zeros_like(noise)buffer[0] = (1 - p) * noise[0]for i in range(1, N): buffer[i] = (1-p)*noise[i] + p*buffer[i-1]noise = buffer
Необходимо всегда создавать еще один массив ради избежание ошибок. Может, здесь его можно было и не использовать, но в следующей фильтре без него не обойтись.
2) Pick-position comb filter (Гребенчатый фильтр) $inline$H_ (z)$inline$.
$$display$$H_ (z)=1-z^{-int(*N+1/2)},(0,1).$$display$$
Соответствующая формула:
$$display$$y(n)=x(n)-x(n-int(*N+1/2)).$$display$$
Код:
pick = int(beta*N+1/2)if pick == 0: pick = N # То есть фильтр не будет действоватьbuffer = np.zeros_like(noise)for i in range(N): if i-pick < 0: buffer[i] = noise[i] else: buffer[i] = noise[i]-noise[i-pick]noise = buffer
В первом абзаце на странице 13 этого документа написано следующее (не дословно, но с сохранением смысла): коэффициент имитирует расположение щипка струны. Если $inline$=1/2$inline$, то это значит, что щипок произвели на середине струны. Если $inline$=1/10$inline$ щипок произвели на одной десятой части струны от моста.
Часть II. Функции, относящиеся к основной части алгоритма
Тут есть ловушка, которую нам придется обойти. Вот, например, String-dampling filter $inline$H_d (z)$inline$ записывается так: $inline$H_d (z)=(1-S)+Sz^{-1}$inline$. Но по рисунку видно, что он берет значение от туда, куда и отдает. То есть получается, что входной и выходной сигналы для этого фильтра это одно и то же. Это означает, что каждый фильтр нельзя применить отдельно, как в прошлой части, надо все фильтры применять одновременно. Это можно сделать, например, найдя произведение каждого фильтра. Но этот подход не рационален: при добавлении или изменении фильтра, придется все снова умножать. Это сделать возможно, но в этом нет смысла. Хотелось бы в один клик менять фильтр, а не умножать все снова и снова.
Так как выходной сигнал от фильтра считается входным для другого фильтра, то я предлагаю написать каждый фильтр отдельной функцией, вызывающей внутри себя функцию прошлого фильтра.
Думаю, на примере кода будет понятно, что я имею в виду.
1) Delay Line filter $inline$z^{-N}. $inline$
$$display$$H(z)=z^{-N}.$$display$$
Соответствующая формула:
$$display$$y(n)=x(n-N).$$display$$
Код:
# Неявно использутся тот факт, что на конце массива samples значение 0.# То есть при n-N<0 значение будет 0, как и должно быть.def DelayLine(n): return samples[n-N]
2) String-dampling filter $inline$H_d (z)$inline$.
$$display$$H_d (z)=(1-S)+Sz^{-1},S[0,1]. $$display$$
В оригинальном алгоритме $inline$S=0.5.$inline$
Соответствующая формула:
$$display$$y(n)=(1-S)x(n)+Sx(n-1).$$display$$
Код:
# String-dampling filter.# H(z) = 0.996*((1-S)+S*z^(-1)). В оригинальном алгоритме S = 0.5. S [0, 1]# y(n)=0.996*((1-S)*x(n)+S*x(n-1))def StringDampling_filter(n): return 0.996*((1-S)*DelayLine(n)+S*DelayLine(n-1))
В данном случае этот фильтр является One Zero String-dampling filter. Существуют и другие варианты, про них можно почитать здесь.
3) String-stiffness allpass filter $inline$H_s (z)$inline$.
Сколько бы я не искал, увы, я не смог найти чего-то конкретного. Здесь этот фильтр написан в общем виде. Но это ничего не дает, так как самая сложная часть это найти подходящие коэффициенты. Еще что-то есть в этом документе на 14 странице, но мне не хватает математической базы, чтобы понять, что там происходит и как это использовать. Если кто-то сможет дайте знать.
4) First-order string-tuning allpass filter $inline$H_ (z)$inline$.
Страница 6, слева внизу в этом документе:
$$display$$H_ (z)=\frac{C+z^{-1}}{1+Cz^{-1}},C(-1,1).$$display$$
Соответствующая формула:
$$display$$y(n)=Cx(n)+x(n-1)-Cy(n-1).$$display$$
Код:
# First-order string-tuning allpass filter# H(z) = (C+z^(-1))/(1+C*z^(-1)). C (-1, 1)# y(n) = C*x(n)+x(n-1)-C*y(n-1)def FirstOrder_stringTuning_allpass_filter(n): # Тут следовало использовать буфер и хранить прошлые значение, но это не нужно, так как # неявно используется тот факт, что прошлое значение уже сохранено и записано в массив samples. return C*(StringDampling_filter(n)-samples[n-1])+StringDampling_filter(n-1)
Нужно помнить, что если вы после этого фильтра добавите еще фильтры, то придется хранить прошлое значение, ибо оно уже не будет хранится в массиве samples.
Так как длина начального шума выражается целым числом, мы выбрасываем дробную часть при подсчете. Это вызывает ошибки и неточности. К примеру, если частота дискретизации равна 44100, а длина шума 133 и 134, то соответствующие значения частоты сигнала равны 331,57 Гц и 329,10 Гц. А частота ноты ми первой октавы (первая открытая струна) 329.63 Гц. Тут разница в десятых, но, например, для 15 лада, разница уже может быть в несколько Гц. Для уменьшения этой ошибки и существует этот фильтр. Его можно не использовать, если частота дискретизации большая (реально большая: несколько сотен тысяц Гц, а то больше) или основная частота мала, как, например, для басовых струн.
Cуществуют и другие вариации, прочитать про них можно все там же.
5) Используем наши функции.
def Modeling(n): return FirstOrder_stringTuning_allpass_filter(n) for i in range(N, len(samples)): samples[i] = Modeling(i)
Часть III. Dynamic Level Lowpass Filter $inline$H_L (z).$inline$
$inline$\check{}=\frac{T}{2}=2f\frac{T}{2}=\frac{f}{F_s}$inline$, где $inline$f$inline$ основная частота, $inline$F_s $inline$ частота дискретизации.
Сначала мы находим массив $inline$y$inline$ следующей формулой:
$$display$$H(z)=\frac{\check{}}{1+\check{}}\frac{1+z^{-1}}{1-\frac{1-\check{}}{1+\check{}} z^{-1}}$$display$$
Соответствующая формула:
$$display$$y(n)=\frac{\check{}}{1+\check{}}(x(n)+x(n-1))+\frac{1-\check{}}{1+\check{}}y(n-1)$$display$$
После применяем следующую формулу:
$$display$$x(n)=L^{\frac{4}{3}}x(n)+(1-L)y(n),L(0,1)$$display$$
Код:
# Dynamic-level lowpass filter. L (0, 1/3)w_tilde = np.pi*frequency/sample_ratebuffer = np.zeros_like(samples)buffer[0] = w_tilde/(1+w_tilde)*samples[0]for i in range(1, len(samples)): buffer[i] = w_tilde/(1+w_tilde)*(samples[i]+samples[i-1])+(1-w_tilde)/(1+w_tilde)*buffer[i-1]samples = (L**(4/3)*samples)+(1.0-L)*buffer
Параметр L влияет на значение уменьшение громкости. При его значениях равных 0.001, 0.01, 0.1, 0.32 громкость сигнала уменьшается на 60, 40, 20 и 10 дБ соответственно.
Оформим все как функцию. Собственно, вот и весь код.
import numpy as npimport scipy.io.wavfile as wave def GuitarString(frequency, duration=1., sample_rate=44100, p=0.9, beta=0.1, S=0.5, C=0.1, L=0.1, toType=False): N = int(sample_rate/frequency) # Длина шума в сэмплах noise = np.random.uniform(-1, 1, N) # Создаем шум # Pick-direction lowpass filter (Фильтр низких частот). # H(z) = (1-p)/(1-p*z^(-1)). p [0, 1) # y(n) = (1-p)*x(n)+p*y(n-1) buffer = np.zeros_like(noise) buffer[0] = (1 - p) * noise[0] for i in range(1, N): buffer[i] = (1-p)*noise[i] + p*buffer[i-1] noise = buffer # Pick-position comb filter (Гребенчатый фильтр). # H(z) = 1-z^(-int(beta*N+1/2)). beta (0, 1) # y(n) = x(n)-x(n-int(beta*N+1/2)) pick = int(beta*N+1/2) if pick == 0: pick = N # То есть фильтр не будет действовать buffer = np.zeros_like(noise) for i in range(N): if i-pick < 0: buffer[i] = noise[i] else: buffer[i] = noise[i]-noise[i-pick] noise = buffer # Добавляем шум в начале. samples = np.zeros(int(sample_rate*duration)) for i in range(N): samples[i] = noise[i] # Неявно использутся тот факт, что на конце массива samples значение 0. # То есть при n-N<0 значение будет 0, как и должно быть. def DelayLine(n): return samples[n-N] # String-dampling filter. # H(z) = 0.996*((1-S)+S*z^(-1)). В оригинальном алгоритме S = 0.5. S [0, 1] # y(n)=0.996*((1-S)*x(n)+S*x(n-1)) def StringDampling_filter(n): return 0.996*((1-S)*DelayLine(n)+S*DelayLine(n-1)) # First-order string-tuning allpass filter # H(z) = (C+z^(-1))/(1+C*z^(-1)). C (-1, 1) # y(n) = C*x(n)+x(n-1)-C*y(n-1) def FirstOrder_stringTuning_allpass_filter(n): # Тут следовало использовать буфер и хранить прошлые значение, но это не нужно, так как # неявно используется тот факт, что прошлое значение уже сохранено и записано в массив samples. return C*(StringDampling_filter(n)-samples[n-1])+StringDampling_filter(n-1) def Modeling(n): return FirstOrder_stringTuning_allpass_filter(n) for i in range(N, len(samples)): samples[i] = Modeling(i) # Dynamic-level lowpass filter. L (0, 1/3) w_tilde = np.pi*frequency/sample_rate buffer = np.zeros_like(samples) buffer[0] = w_tilde/(1+w_tilde)*samples[0] for i in range(1, len(samples)): buffer[i] = w_tilde/(1+w_tilde)*(samples[i]+samples[i-1])+(1-w_tilde)/(1+w_tilde)*buffer[i-1] samples = (L**(4/3)*samples)+(1.0-L)*buffer if toType: samples = samples/np.max(np.abs(samples)) # Нормируем от -1 до 1 return np.int16(samples*32767) # Переводим в тип данных int16 else: return samples frequency = 82.51sound = GuitarString(frequency, duration=4, toType=True)wave.write("SoundGuitarString.wav", 44100, sound)
Открытая шестая струна (82.41 Гц) звучит так:
А открытая первая струна (329.63 Гц) звучит так:
Первая струна звучит, мягко говоря, не очень. Больше похоже на колокол, чем на струну. Я очень долго пытался понять, что не так в алгоритме. Думал, что дело в неиспользованном фильтре. Спустя дни экспериментов, я понял, что нужно увеличить частоту дискретизации хотя бы до 100000:
Звучит лучше, не правда ли?
Про дополнения, такие как игра глиссандо или симуляция симпатической струны, можно почитать в этом документе (с. 11-12).
Вот вам бой:
Последовательность аккордов: C G Am F. Бой: шестерка. Задержка между двумя последовательными щипками струны равна 0.015 секунд; задержка между двумя последовательными ударами в бою равна 0.41 секунда; сама задержка в бою равна 0.82 секунды. В алгоритме изменено значение L на 0.2.
Напоследок, вот вам наипростейший перебор:
Спасибо за прочтение статьи. Удачи!