Моделирование звука гитарных нот с помощью алгоритма Карплуса-Стронга на python

Знакомьтесь, эталонная нота ля первой октавы (440 Гц):


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

Но нас интересует другой вопрос: как этот уникальный тембр смоделировать на компьютере?

Примечание
В этой статье не будет разбираться почему это работает. Будут лишь ответы на вопросы: что это и как это работает?


Стандартный алгоритм Карплуса-Стронга


image

Иллюстрация взята с этого сайта.

Суть алгоритма в следующем:

1) Создаем массив размера N из случайных чисел (N напрямую связана с основной частотой звука).

2) Добавляем к концу этого массива значение, посчитанное по следующей формуле:

$y(n)=\frac{y(n-N)+y(n-N-1)}{2},$


где $y$ – наш массив.

3) Выполняем пункт 2 необходимое количество раз.

Приступим к написанию кода:

1) Импортируем необходимые библиотеки.

import numpy as np
import 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 np
import 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.41
sound = GuitarString(frequency, duration=4, toType=True)
wave.write("SoundGuitarString.wav", 44100, sound)

9) Запустим и получим:


Для того, чтобы струна звучала лучше, слегка улучшим формулу:

$y(n)=0.996\frac{y(n-N)+y(n-N-1)}{2}$



Открытая шестая струна (82.41 Гц) звучит так:


Открытая первая струна (329.63 Гц) звучит так:


Звучит неплохо, не правда ли?

Можно бесконечно подбирать этот коэффициент и найти среднее между красивым звучанием и длительностью, но лучше сразу перейти к Расширенному алгоритму Карплуса-Стронга.

Немного о Z-преобразовании


Примечание
Этот раздел существует лишь из-за того, что все функции записываются в виде Z-преобразования. Для того, чтобы вы могли смогли воспользоваться формулой или формулами, которые не описаны здесь (а такие существуют), и тем самым смогли улучшить алгоритм, надо понять, как это Z-преобразование использовать. В этом разделе не будут ответы на вопросы: что это, почему и как работает?

Пусть $x$ – массив входных значений, а $y$ – массив выходных значений. Каждый элемент массива y выражается следующей формулой:

$y(n)=x(n)+x(n-1).$



Если индекс выходит за пределы массива, то значение равно 0. То есть $x(0-1)=0$. (Посмотрите прошлый код, там это неявно использовалось).

Эту формулу можно записать в соответствующем Z-преобразовании:

$H(z)=1+z^{-1}.$



Если формула такая:

$y(n)=x(n)+x(n-1)-y(n-1).$



То есть каждый элемент входного массива зависит от прошлого элемента этого же массива (кроме нулевого элемента, конечно). То соответствующее Z-преобразование выглядит так:

$H(z)=\frac{1+z^{-1}}{1+z^{-1}}.$


Обратный процесс: из Z-преобразования получить формулу для каждого элемента. Например,

$H(z)=\frac{1+z^{-1}}{1-z^{-1}}.$


$H(z)=\frac{Y(z)}{X(z)} =\frac{1+z^{-1}}{1-z^{-1}}. $


$Y(z)*(1-z^{-1} )=X(z)*(1+z^{-1} ). $


$Y(z)*1-Y(z)*z^{-1}= X(z)*1+X(z)*z^{-1}.$


$y(n)-y(n-1)=x(n)+x(n-1). $


$y(n)=x(n)+x(n-1)+y(n-1).$


Если кто-то не понял, то формула такая: $Y(z)*α*z^{-k}=α*y(n-k)$, где $α$ – любое действительное число.

Если надо умножить два Z-преобразования друг на друга, то $z^{-a}*z^{-b}=z^{-a-b}.$

Расширенный алгоритм Карплуса-Стронга


image
Иллюстрация взята с этого сайта.

Вот быстрое описание каждой функции.

Часть I. Функции, преобразующие начальный шум


1) Pick-direction lowpass filter (Фильтр низких частот) $H_p (z)$.

$H_p (z)=\frac{1-p}{1-pz^{-1}},p ∈ [0,1).$


Соответствующая формула:

$y(n)=(1-p)x(n)+py(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

Необходимо всегда создавать еще один массив ради избежание ошибок. Может, здесь его можно было и не использовать, но в следующей фильтре без него не обойтись.

2) Pick-position comb filter (Гребенчатый фильтр) $H_β (z)$.

$$display$$H_β (z)=1-z^{-int(β*N+1/2)},β∈(0,1).$$display$$


Соответствующая формула:

$y(n)=x(n)-x(n-int(β*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

В первом абзаце на странице 13 этого документа написано следующее (не дословно, но с сохранением смысла): коэффициент β имитирует расположение щипка струны. Если $β=1/2$, то это значит, что щипок произвели на середине струны. Если $β=1/10$ — щипок произвели на одной десятой части струны от моста.

Часть II. Функции, относящиеся к основной части алгоритма


Тут есть ловушка, которую нам придется обойти. Вот, например, String-dampling filter $H_d (z)$ записывается так: $H_d (z)=(1-S)+Sz^{-1}$. Но по рисунку видно, что он берет значение от туда, куда и отдает. То есть получается, что входной и выходной сигналы для этого фильтра это одно и то же. Это означает, что каждый фильтр нельзя применить отдельно, как в прошлой части, надо все фильтры применять одновременно. Это можно сделать, например, найдя произведение каждого фильтра. Но этот подход не рационален: при добавлении или изменении фильтра, придется все снова умножать. Это сделать возможно, но в этом нет смысла. Хотелось бы в один клик менять фильтр, а не умножать все снова и снова.
Так как выходной сигнал от фильтра считается входным для другого фильтра, то я предлагаю написать каждый фильтр отдельной функцией, вызывающей внутри себя функцию прошлого фильтра.
Думаю, на примере кода будет понятно, что я имею в виду.
1) Delay Line filter $z^{-N}. $

$H(z)=z^{-N}.$


Соответствующая формула:

$y(n)=x(n-N).$


Код:

# Неявно использутся тот факт, что на конце массива samples значение 0.
# То есть при n-N<0 значение будет 0, как и должно быть.
def DelayLine(n):
    return samples[n-N]


2) String-dampling filter $H_d (z)$.

$H_d (z)=(1-S)+Sz^{-1},S∈[0,1]. $


В оригинальном алгоритме $S=0.5.$
Соответствующая формула:

$y(n)=(1-S)x(n)+Sx(n-1).$


Код:

# 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 $H_s (z)$.
Сколько бы я не искал, увы, я не смог найти чего-то конкретного. Здесь этот фильтр написан в общем виде. Но это ничего не дает, так как самая сложная часть – это найти подходящие коэффициенты. Еще что-то есть в этом документе на 14 странице, но мне не хватает математической базы, чтобы понять, что там происходит и как это использовать. Если кто-то сможет – дайте знать.

4) First-order string-tuning allpass filter $H_ρ (z)$.
Страница 6, слева внизу в этом документе:

$$display$$H_ρ (z)=\frac{C+z^{-1}}{1+Cz^{-1}},C∈(-1,1).$$display$$


Соответствующая формула:

$y(n)=Cx(n)+x(n-1)-Cy(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)

Нужно помнить, что если вы после этого фильтра добавите еще фильтры, то придется хранить прошлое значение, ибо оно уже не будет хранится в массиве 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 $H_L (z).$


$inline$\check{ω}=ω\frac{T}{2}=2πf\frac{T}{2}=π\frac{f}{F_s}$inline$, где $f$ – основная частота, $F_s $– частота дискретизации.
Сначала мы находим массив $y$ следующей формулой:

$$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$$


После применяем следующую формулу:

$x(n)=L^{\frac{4}{3}}x(n)+(1-L)y(n),L∈(0,1)$


Код:

# 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

Параметр L влияет на значение уменьшение громкости. При его значениях равных 0.001, 0.01, 0.1, 0.32 громкость сигнала уменьшается на 60, 40, 20 и 10 дБ соответственно.

Оформим все как функцию. Собственно, вот и весь код.

import numpy as np
import 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.51
sound = 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.205 секунда; сама задержка в бою равна 0.41 секунды. В алгоритме изменено значение L на 0.2.

Спасибо за прочтение статьи. Удачи!
AdBlock похитил этот баннер, но баннеры не зубы — отрастут

Подробнее
Реклама

Комментарии 7

    0
    Класс! А не пробовали сравнивать сигнал с реальной гитары и синтезированный?
    Может можно подобрать коэффициенты фильтров (например, генетическим алгоритмом), чтобы добиться наибольшей «похожести»?
      0

      Ваш коммент натолкнул на мысль, что если в ИИ на вход подавать синтезированный сингал а в качестве правильного ответа вживую сыгранный.
      Сможет ли она потом сама его преобразовывать?

        0
        Скорее всего да. Только по скорости вычисления будет хуже.
        Можно пойти еще дальше. Подавать на вход шум и номер ноты, а на выходе получать нужный сигнал. Возможно нейросеть подберет «правильные» фильтры и задержки.
          0
          Этак вы и до WaveNet дойдёте.
        0
        Идея с генетическим алгоритмом мне нравится!
        В качестве сравнения сигналов предлагаю сравнивать их на фундаментальном уровне: сравнивать значение гармоник с наибольшими амплитудами, использовав быстрое преобразование Фурье. А шум и погрешность можно исключить, например, находя максимумы, только если амплитуда гармоники выше, чем средняя амплитуда всех гармоник. Вот в этой статье хорошо видны максимумы гармоник.
          +1
          Думаю гармоники не совсем то что нужно.
          Когда-то давно я тоже пробовал сигналы с гитары обрабатывать, но методом Прони.
          Сам затухающий синус с гармониками воспроизводится хорошо. Для этого не так много гармоник требуется. Но если мы отбрасываем какие-либо гармоники, то потеряем начальный удар о струну.
          На мой взгляд именно его труднее всего воспроизвести.
        +5

        Все-таки какой-то грязновато-звенящий бой получился. Меня в свое время впечатлила реализация Андре Мишеля (это на флеше было) — там сразу несколько параметров типа приглушения и т.д.


        Его демку позже на джаваскрипте переписали, послушайте: http://amid.fish/javascript-karplus-strong

        Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

        Самое читаемое