pythonのKarplus-Strongアルゴリズムを使用したギターノートのサウンドのモデリング

最初のオクターブ（440 Hz）のリファレンスノートAを満たします。





痛そうですね。同じ音が異なる楽器で異なって聞こえるという事実について他に何を言うべきか。なぜそうなのですか？それはすべて、各楽器に固有の音色を作り出す追加の高調波の存在に関するものです。



しかし、私たちは別の質問に興味があります：コンピュータでこのユニークな音色をシミュレートする方法は？

標準のKarplus-強力なアルゴリズム

このサイトから取られたイラスト。



アルゴリズムの本質は次のとおりです



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

開いた6番目の文字列（82.41 Hz）は次のように聞こえます。





開いた最初の文字列（329.63 Hz）は次のように聞こえます。





いいですね。



この係数を際限なく選択して、美しいサウンドと持続時間の平均を見つけることができますが、AdvancedKarplus-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)$どこ $$$α$-任意の実数。



2つのZ変換を互いに乗算する必要がある場合は、$$$z^{-a}*z^{-b}=z^{-a-b}.$

拡張Karplus-強力なアルゴリズム

このサイト から取られたイラスト。



ここだ各機能の簡単な説明が。

パートI.初期ノイズを変換する関数

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$ -橋からの弦の10分の1を引っ張った。

パートII。アルゴリズムの主要部分に関連する関数

ここには回避しなければならない罠があります。たとえば、文字列ダンピングフィルタ$$$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))

この場合、このフィルターはOne ZeroString-damplingフィルターです。他のオプションがあります、あなたはここでそれらについて読むことができます。



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.57Hzと329.10Hzです。そして、最初のオクターブ（最初の開いたストリング）のEノートの周波数は329.63Hzです。ここでは、差は10分の1ですが、たとえば、15フレットの場合、差はすでに数Hzになっている可能性があります。このエラーを減らすために、このフィルターが存在します。サンプリング周波数が高い場合（実際には数十万Hz以上）、または低音ストリングのように基本周波数が低い場合は省略できます。

他のバリエーションがありますが、あなたはすべてのそれらについて読むことができますが。



5）機能を使用します。

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

パートIII。ダイナミックレベルローパスフィルター $$$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、および10dB減少します。



すべてを関数として設計しましょう。実際には、それがすべてのコードです。

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)

開いた6番目の文字列（82.41 Hz）は次のように聞こえます。





そして、開いた最初の文字列（329.63 Hz）は次のように聞こえます。





最初の弦は、穏やかに言えば、あまり良く聞こえません。ひもというよりは鐘のようなものです。私は非常に長い間、アルゴリズムの何が問題なのかを理解しようとしてきました。未使用のフィルターだと思った。数日間の実験の後、サンプルレートを少なくとも100,000に増やす必要があることに気付きました。





良く聞こえますね。このドキュメント



では、グリサンドの再生や同情的な文字列のシミュレーションなどのアドオンを読むことができます（11〜12ページ）。 これがあなたのための戦いです：









コード進行：CG＃Am F.ストライク：6。ストリングを2回連続して引き抜く間の遅延は0.015秒です。戦闘での2つの連続したヒット間の遅延は0.205秒です。戦闘の遅延自体は0.41秒です。アルゴリズムはLの値を0.2に変更しました。



記事を読んでいただきありがとうございます。幸運を！