Modelando o som de notas de guitarra usando o algoritmo Karplus-Strong em python

Conheça a nota de referência A da primeira oitava (440 Hz):





Parece doloroso, não é? O que mais dizer sobre o fato de que a mesma nota soa de maneira diferente em diferentes instrumentos musicais. Por que é tão? É tudo sobre a presença de harmônicos adicionais que criam um timbre único para cada instrumento.



Mas estamos interessados ​​em outra questão: como simular esse timbre único em um computador?

Algoritmo Karplus-Strong padrão

Ilustração retirada deste site .



A essência do algoritmo é a seguinte:



1) Crie uma matriz de tamanho N a partir de números aleatórios (N está diretamente relacionado à frequência de som fundamental).



2) Adicione ao final desta matriz o valor calculado pela seguinte fórmula:

$$

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

Onde $$$y$É a nossa matriz.



3) Realizamos o ponto 2 o número necessário de vezes.



Vamos começar a escrever o código:



1) Importe as bibliotecas necessárias.

import numpy as np
import scipy.io.wavfile as wave

2) Inicializamos as variáveis.

frequency = 82.41     #     
duration = 1          #    
sample_rate = 44100   #  

3) Crie ruído.

#  ,  frequency, ,        frequency .
#      sample_rate/length .
#  length = sample_rate/frequency.
noise = np.random.uniform(-1, 1, int(sample_rate/frequency))   

4) Crie uma matriz para armazenar os valores e adicionar ruído no início.

samples = np.zeros(int(sample_rate*duration))
for i in range(len(noise)):
    samples[i] = noise[i]

5) Usamos a fórmula.

for i in range(len(noise), len(samples)):
    #   i   ,      .
    #  ,  i   ,       .
    samples[i] = (samples[i-len(noise)]+samples[i-len(noise)-1])/2

6) Normalizamos e traduzimos para o tipo de dados desejado.

samples = samples / np.max(np.abs(samples))  
samples = np.int16(samples * 32767)     

7) Salvar em arquivo.

wave.write("SoundGuitarString.wav", 44100, samples)

8) Vamos projetar tudo como uma função. Na verdade, esse é todo o código.

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) Vamos correr e obter:





Para fazer a corda soar melhor, vamos melhorar um pouco a fórmula:

$$

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

Uma sexta corda aberta (82,41 Hz) soa assim:





A primeira corda aberta (329,63 Hz) soa assim:





Parece bom, não é?



Você pode selecionar esse coeficiente indefinidamente e encontrar a média entre o som bonito e a duração, mas é melhor ir direto para o algoritmo Karplus-Strong avançado.

Um pouco sobre a transformada Z

Deixe ser $$$x$ É uma matriz de valores de entrada e $$$y$- uma matriz de valores de saída. Cada elemento em y é expresso pela seguinte fórmula:

$$

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

Se o índice estiver fora da matriz, o valor é 0. Isso é $$$x(0-1)=0$... (Veja o código anterior, lá foi usado implicitamente).



Esta fórmula pode ser escrita na transformação Z correspondente:

$$

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

Se a fórmula for assim:

$$

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

Ou seja, cada elemento da matriz de entrada depende do elemento anterior da mesma matriz (exceto para o elemento zero, é claro). Então, a transformação Z correspondente se parece com isto:

$$

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

Processo reverso: obtenha a fórmula para cada elemento da transformada Z. Por exemplo,

$$

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

Se alguém não entende, a fórmula é: $$$Y(z)*α*z^{-k}=α*y(n-k)$Onde $$$α$- qualquer número real.



Se você precisar multiplicar duas transformações Z uma pela outra, então$$$z^{-a}*z^{-b}=z^{-a-b}.$

Algoritmo Karplus-Strong estendido

Ilustração retirada deste site.



Aqui está um rápido resumo de cada recurso.

Parte I. Funções que transformam o ruído inicial

1) Filtro passa- baixo de direção de seleção ( filtro passa- baixo)$$$H_p (z)$...

$$

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

Fórmula correspondente:

$$

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

O código:

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

Você deve sempre criar outro array para evitar erros. Talvez não pudesse ter sido usado aqui, mas no próximo filtro você não pode fazer sem ele.



2) Filtro pente de posição de seleção (filtro pente)$$$H_β (z)$...

$$

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

Fórmula correspondente:

$$

$$y(n)=x(n)-x(n-int(β*N+1/2)).$$

O código:

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

No primeiro parágrafo da página 13 deste documento, está escrito o seguinte (não literalmente, mas com a preservação do significado): o coeficiente β imita a posição da corda puxada. E se$$$β=1/2$, significa que a dedilhada foi feita no meio da corda. E se$$$β=1/10$ - a arrancada foi feita em um décimo da corda da ponte.

Parte II. Funções relacionadas à parte principal do algoritmo

Há uma armadilha aqui que temos que contornar. Por exemplo, filtro de amortecimento de corda$$$H_d (z)$ escrito assim: $$$H_d (z)=(1-S)+Sz^{-1}$... Mas a figura mostra que ele tira sentido de onde o dá. Ou seja, verifica-se que os sinais de entrada e saída para este filtro são um e o mesmo. Isso significa que cada filtro não pode ser aplicado separadamente, como na seção anterior, todos os filtros devem ser aplicados simultaneamente. Isso pode ser feito, por exemplo, encontrando o produto de cada filtro. Mas essa abordagem não é racional: ao adicionar ou alterar um filtro, você terá que multiplicar tudo novamente. É possível fazer isso, mas não faz sentido. Eu gostaria de mudar o filtro em um clique e não multiplicar tudo repetidamente.

Visto que o sinal de saída do filtro é considerado a entrada de outro filtro, proponho escrever cada filtro como uma função separada que chama a função do filtro anterior dentro de si.

Acho que o código de exemplo deixará claro o que quero dizer.

1) Filtro de linha de atraso $$$z^{-N}.$

$$

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

Fórmula correspondente:

$$

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

O código:

#    ,     samples  0.
#    n-N<0   0,    .
def DelayLine(n):
    return samples[n-N]

2) Filtro de amortecimento de corda $$$H_d (z)$...

$$

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

No algoritmo original $$$S=0.5.$

Fórmula correspondente:

$$

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

O código:

# 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))

Nesse caso, esse filtro é o filtro One Zero String-dampling. Existem outras opções, você pode ler sobre elas aqui .



3) Filtro passa-tudo de rigidez de corda$$$H_s (z)$...

Por mais que procurasse, infelizmente, não consegui encontrar nada específico. Aqui, este filtro é escrito em termos gerais. Mas isso não funciona, pois a parte mais difícil é encontrar as chances certas. Há algo mais neste documento na página 14, mas não tenho conhecimento matemático suficiente para entender o que está acontecendo lá e como usá-lo. Se alguém puder, me avise.



4) Filtro passa-tudo de sintonia de primeira ordem$$$H_ρ (z)$...

Página 6, canto inferior esquerdo deste documento:

$$

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

Fórmula correspondente:

$$

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

O código:

# 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)

Deve ser lembrado que se você adicionar mais filtros após este filtro, você terá que armazenar o valor anterior, porque ele não será mais armazenado no array samples.

Como o comprimento do ruído inicial é um inteiro, descartamos a parte fracionária ao contar. Isso causa erros e imprecisões. Por exemplo, se a taxa de amostragem for 44100 e o comprimento do ruído for 133 e 134, então as frequências de sinal correspondentes são 331,57 Hz e 329,10 Hz. E a frequência das notas E da primeira oitava (a primeira corda aberta) é 329,63 Hz. Aqui a diferença está em décimos, mas, por exemplo, para o 15º traste, a diferença já pode ser de vários Hz. Este filtro existe para reduzir este erro. Pode ser omitido se a frequência de amostragem for alta (realmente alta: várias centenas de milhares de Hz, ou até mais) ou a frequência fundamental for baixa, como, por exemplo, para cordas graves.

Existem outras variações, você pode ler sobre todas elas lá .



5) Usamos nossas funções.

def Modeling(n):
    return FirstOrder_stringTuning_allpass_filter(n)
 
for i in range(N, len(samples)):
    samples[i] = Modeling(i)

Parte III. Filtro Lowpass de nível dinâmico $$$H_L (z).$

$$$\check{ω}=ω\frac{T}{2}=2πf\frac{T}{2}=π\frac{f}{F_s}$Onde $$$f$ - frequência fundamental, $$$F_s$- frequência de amostragem.

Primeiro encontramos a matriz$$$y$ com a seguinte fórmula:

$$

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

Fórmula correspondente:

$$

$$y(n)=\frac{\check{ω}}{1+\check{ω}}(x(n)+x(n-1))+\frac{1-\check{ω}}{1+\check{ω}}y(n-1)$$

Em seguida, aplicamos a seguinte fórmula:

$$

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

O código:

# 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

O parâmetro L afeta o valor de diminuição do volume. Com seus valores iguais a 0,001, 0,01, 0,1, 0,32, o volume do sinal diminui em 60, 40, 20 e 10 dB, respectivamente.



Vamos projetar tudo como uma função. Na verdade, esse é todo o código.

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)

Uma sexta corda aberta (82,41 Hz) soa assim:





E a primeira corda aberta (329,63 Hz) soa assim:





A primeira corda não soa muito bem, para dizer o mínimo. Mais como um sino do que uma corda. Por muito tempo tentei descobrir o que havia de errado no algoritmo. Pensei que fosse um filtro não utilizado. Após dias de experiências, percebi que precisava aumentar a taxa de amostragem para pelo menos 100.000:





Parece melhor, não é?



Complementos como tocar glissando ou simular uma string simpática podem ser lidos neste documento (pp. 11-12).



Aqui está uma luta:





Sequência de acordes: CG # Am F. Strike: Seis. O atraso entre duas arrancadas consecutivas da corda é de 0,015 segundos; o atraso entre dois acertos consecutivos em uma batalha é de 0,205 segundos; o próprio atraso na batalha é de 0,41 segundos. O algoritmo mudou o valor de L para 0,2.



Obrigado por ler o artigo. Boa sorte!