نمذجة صوت نغمات الجيتار باستخدام خوارزمية Karplus-Strong في Python

تعرف على الملاحظة المرجعية أ لأول أوكتاف (440 هرتز):





تبدو مؤلمة ، أليس كذلك؟ ماذا نقول عن حقيقة أن نفس النوتة تبدو مختلفة على الآلات الموسيقية المختلفة. لماذا هو كذلك؟ الأمر كله يتعلق بوجود التوافقيات الإضافية التي تخلق جرسًا فريدًا لكل أداة.



لكننا مهتمون بسؤال آخر: كيف نحاكي هذا الجرس الفريد على جهاز كمبيوتر؟

خوارزمية 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 هرتز) تبدو كالتالي:





تبدو السلسلة الأولى المفتوحة (329.63 هرتز) كما يلي:





تبدو جيدة ، أليس كذلك؟



يمكنك تحديد هذا المعامل إلى ما لا نهاية والعثور على المتوسط ​​بين الصوت الجميل والمدة ، ولكن من الأفضل الانتقال مباشرة إلى خوارزمية 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-transform. على سبيل المثال،

$$

$$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) مرشح مشط Pick-position (مرشح مشط)$$$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$ - تم عمل النتف على عُشر الخيط من الجسر.

الجزء الثاني. الوظائف المتعلقة بالجزء الرئيسي من الخوارزمية

هناك فخ هنا علينا أن نتغلب عليه. على سبيل المثال ، مرشح String-dampling$$$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 Zero String-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.57 هرتز و 329.10 هرتز. وتردد النغمات E لأول أوكتاف (السلسلة المفتوحة الأولى) هو 329.63 هرتز. الفرق هنا هو بالأعشار ، ولكن ، على سبيل المثال ، بالنسبة للحنق الخامس عشر ، قد يكون الفرق بالفعل عدة هرتز. لتقليل هذا الخطأ ، يوجد عامل التصفية هذا. لا يجوز استخدامه إذا كان تردد أخذ العينات مرتفعًا (مرتفع حقًا: عدة مئات من آلاف هرتز ، أو حتى أكثر) أو إذا كان التردد الأساسي منخفضًا ، على سبيل المثال ، لسلاسل الجهير.

هناك اختلافات أخرى، يمكنك أن تقرأ عن كل منهم هناك .



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 ديسيبل على التوالي.



دعونا نصمم كل شيء كوظيفة. في الواقع ، هذا كل ما في الكود.

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 على الأقل:





يبدو أفضل ، أليس كذلك؟



يمكن قراءة الوظائف الإضافية مثل عزف glissando أو محاكاة سلسلة متعاطفة في هذا المستند (الصفحات 11-12).



هنا قتال:





تسلسل الوتر: CG # Am F. Strike: Six. التأخير بين نتف السلسلة مرتين متتاليتين هو 0.015 ثانية ؛ التأخير بين ضربتين متتاليتين في معركة هو 0.205 ثانية ؛ التأخير نفسه في المعركة 0.41 ثانية. لقد غيرت الخوارزمية قيمة L إلى 0.2.



شكرا لقرائتك المجلة. حظا سعيدا!