使用python中的Karplus-Strong算法为吉他音符建模

符合第一个八度音阶（440 Hz）的参考音符A：





听起来很痛苦，不是吗？关于同一音符在不同乐器上听起来不同的事实还有什么要说的。为什么会这样呢？一切都与附加谐波的存在有关，这些附加谐波会为每台乐器创建一个独特的音色。



但是我们对另一个问题感兴趣：如何在计算机上模拟这种独特的音色？

标准Karplus-Strong算法

插图取自本网站。



该算法的本质如下：



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 Hz）听起来像这样：





打开的第一弦（329.63 Hz）听起来像这样：





听起来不错，不是吗？



您可以无休止地选择该系数，并找到优美声音和持续时间之间的平均值，但是最好直接使用Advanced Karplus-Strong算法。

关于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}.$

扩展Karplus-Strong算法

该图取自本网站。



这是每个功能的快速摘要。

第一部分：转换初始噪声的函数

1）拾取方向低通滤波器（低通滤波器）$$$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）取位梳状过滤器（梳状过滤器）$$$H_β (z)$...

$$

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

对应公式：

$$

$$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$ -从琴桥的细绳的十分之一拔出。

第二部分 与算法主要部分相关的功能

这里有一个陷阱，我们必须解决。例如，字符串采样过滤器$$$H_d (z)$ 这样写： $$$H_d (z)=(1-S)+Sz^{-1}$... 但是图片显示，他从他给予的意义中汲取了意义。也就是说，事实证明该滤波器的输入和输出信号是相同的。这意味着不能单独应用每个过滤器，如上一节所述，必须同时应用所有过滤器。例如，这可以通过找到每个过滤器的乘积来完成。但是这种方法并不合理：添加或更改过滤器时，您将不得不再次将所有内容相乘。可以这样做，但这没有任何意义。我想一键更改过滤器，而不是一遍又一遍地乘以所有内容。

由于来自滤波器的输出信号被视为另一个滤波器的输入，因此我建议将每个滤波器写为一个单独的函数，该函数在内部调用前一个滤波器的函数。

我认为示例代码将清楚说明我的意思。

1）延迟线滤波器 $$$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）字符串采样滤波器 $$$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))

在这种情况下，此过滤器是“一零字符串采样”过滤器。还有其他选择，您可以在此处阅读有关它们的信息。



3）弦刚度全通滤波器$$$H_s (z)$...

不管我看上去如何，,，我找不到任何具体的东西。这里的过滤器是用一般术语写的。但这并不可行，因为最困难的部分是找到合适的赔率。该文档在第14页中还有其他内容，但是我没有足够的数学基础来了解那里正在发生的事情以及如何使用它。如果有人可以，请告诉我。



4）一阶调弦全通滤波器$$$H_ρ (z)$... 本文档

左下方第6页：

$$

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

对应公式：

$$

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

应当记住，如果在此过滤器之后添加更多过滤器，则必须存储过去的值，因为它将不再存储在示例数组中。

由于初始噪声的长度是整数，因此在计数时我们将小数部分丢弃。这会导致错误和不准确。例如，如果采样率为44100，噪声长度为133和134，则相应的信号频率为331.57 Hz和329.10 Hz。第一个八度的E音符（第一个打开的弦）的频率为329.63 Hz。在这里，差异为十分之一，但是，例如，对于第15品格，差异可能已经是几赫兹。为减少此错误，此过滤器存在。如果采样频率很高（真的很高：几十万赫兹，甚至更高），或者基频很低，例如低音弦，则可以省略。

还有其他变体，您可以在那里全部阅读。



5）我们使用我们的功能。

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

第三部分 动态电平低通滤波器 $$$H_L (z).$

$$$\check{ω}=ω\frac{T}{2}=2πf\frac{T}{2}=π\frac{f}{F_s}$哪里 $$$f$ -基本频率 $$$F_s$- 采样频率。

首先我们找到数组$$$y$ 具有以下公式：

$$

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

对应公式：

$$

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

然后我们应用以下公式：

$$

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



让我们将所有功能设计为功能。实际上，这就是所有代码。

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 Hz）听起来像这样：





打开的第一个字符串（329.63 Hz）听起来像这样：





轻轻地说，第一弦听起来不太好。更像是铃而不是弦。很长时间以来，我一直试图找出算法中的错误。以为它是一个未使用的过滤器。经过几天的试验，我意识到我需要将采样率提高到至少100,000：





听起来更好，不是吗？



可以在本文档中阅读诸如玩滑音或模拟有同感的弦乐之类的附件（第11-12页）。



这是为您而战：





和弦进行：CG＃AmF。打击：六。两次连续拔弦之间的延迟为0.015秒；战斗中两次连续命中之间的延迟时间为0.205秒；战斗本身的延迟是0.41秒。该算法已将L的值更改为0.2。



感谢您阅读本文。祝好运！