Karplus-Strong 算法简单介绍和实现

• 本文为Coursera数字信号处理课程第一周内容，对相关课程刚兴趣的同学，请参看这里
• 为了有更好的交互性，本文所有代码均上传至Microsoft Azure Notebooks，你可以在上面试听所有输出的音频，具体代码在这里

什么是Karplus-Strong算法

Karplus-Strong算法简单的说就是一种合成声音的算法，它不断的循环重复$M$M个采样点，从而生成一段音频

是的，它非常简单，我们用一段代码（KS_1）来直观的描述它。

import numpy as np
import matplotlib.pyplot as plt
import IPython
%matplotlib inline

def KS_1(x, N):
    '''
    x: 初始序列
    N: 输出序列的最大长度
    '''
    y = x
    while len(y) < N:
        y = np.append(y, x)
        
    y = y[0:N+1]
    
    return y

Fs = 16000 # 设置采样频率，16Khz

当设定了采样频率，那么我们就可以就算用KS算法生成音频的频率了。例如，如果初始序列x的长度为50，采样频率Fs=16000，那么生成音频在每秒内将会重复
16000/50 = 320次x，也就是说生成音频的频率为320Hz，这刚好是钢琴中E4的音。

接下来我们讨论初始序列x的值应该是怎样的，这是KS算法中最coool的地方，因为算法的发明者推荐我们使用高斯白噪声作为初始化序列，虽然我们可以使用任意值的初始化序列，但实验表明，随机的序列效果通常最好

b = np.random.randn(50) # 一个随机的初始化序列
plt.stem(b)

让我们用KS_1来生成一段2s的音频

y = KS_1(b, Fs*2)

plt.figure(figsize=(10,3))
plt.stem(y[0:Fs*2:111])
plt.show()

IPython.display.Audio(y, rate=Fs) # 真的很难听！！

 # 生成另一段音频，但是初始序列更长，意味着输出音频的频率更小
IPython.display.Audio(KS_1(np.random.rand(100), Fs*2), rate=Fs)

很棒！KS算法至少能够运行了。我们可以用下面的图来描述这个信号处理系统

输出信号可以表示为（暂时忽略那个$\alpha$α）
$y[n] = x[n] + y[n-M]$y[n]=x[n]+y[n−M]
并假设x是 finite-support 信号
$x[n] = \begin{cases} 0 &amp; n&lt; 0\\ b_n &amp; 0 \le n &lt; M \\ 0 &amp; n \ge M\\ \end{cases}$x[n]=⎩⎪⎨⎪⎧​0bn​0​n<00≤n<Mn≥M​

让我们上面的式子来实现KS（不是高效的做法）

def KS_2(x, N):
    y = np.zeros(N)
    M = len(x)
    
    for n in range(N):
        y[n] = (x[n] if n < M else 0) + (y[n-M] if n-M>=0 else 0)
        
    return y

IPython.display.Audio(KS_2(np.random.randn(50), Fs*2), rate=Fs)

我们只要添加一个$\alpha$α就能让输出的音频听上去更真实，$\alpha$α设置为比1小那么一丢丢，那么生成的音频听上去就像吉他一样

def KS_3(x, N, alpha=0.99):
    y = np.zeros(N)
    M = len(x)
    
    for n in range(N):
        y[n] = (x[n] if n < M else 0) + alpha*(y[n-M] if n-M>=0 else 0)
        
    return y

y = KS_3(b, Fs*2, 0.99)

plt.figure(figsize=(15,4))
plt.stem(y[0:Fs*2:43]) 
plt.show()

IPython.display.Audio(y, rate=Fs) # 听着不错

KS算法中的还有一个重要的细节。每次循环初始序列都会乘上$\alpha$α，那么我们可以将KS写成：
$y[n] = \alpha^{\lfloor n/M \rfloor}x[n \mod M]$y[n]=α⌊n/M⌋x[nmodM]
也就是说当N固定时（输出音频长度固定），输出音频的减弱速度与$\alpha$α和M有关系。换句话说，如果$\alpha$α固定，M越小（音调越高），那么输出的音频减弱的越快

IPython.display.Audio(KS_3(np.random.rand(50), Fs*2), rate=Fs)

y = KS_3(np.random.randn(50), Fs*2) #衰退速率较慢

plt.figure(figsize=(15,4))
plt.stem(y[0:Fs*2:33])

<Container object of 3 artists>

IPython.display.Audio(KS_3(np.random.rand(10), Fs*2), rate=Fs)

y = KS_3(np.random.randn(10), Fs*2) # 衰退速率快

plt.figure(figsize=(15,4))
plt.stem(y[0:Fs*2:33]) 

<Container object of 3 artists>

这样是不好的，我们需要进行改进。我们希望当$\alpha$α一样时，减弱的速率也是一样的。于是，我们最终版的KS算法来了

def KS(x, N, alpha=0.99, REF_LEN=50):
    
    # REF_LEN 越大衰减越慢（感觉尾音越长），
    # 例如REF_LEN=50，表示当 alpha 固定时，你希望衰减速率同M=50一样

    M = len(x)
    beta = alpha ** (float(M) / REF_LEN) # 表示 \beta^n = \alpha，其中n=REF_LEN/M，求得beta
    
    y = np.zeros(N)
    
    for n in range(N):
        y[n] = (x[n] if n < M else 0) + beta*(y[n-M] if n-M >= 0 else 0)
        
    return y

y = KS(np.random.rand(50), Fs*2, REF_LEN=50)
IPython.display.Audio(y, rate=Fs)

y = KS(np.random.rand(10), Fs*2, REF_LEN=50)
IPython.display.Audio(y, rate=Fs)

上面的KS算法效率较低，用上numpy，可以大大提高效率

def KS_with_np(x, N, alpha=0.99, REF_LEN=50):
    M = len(x)
    beta = alpha ** (M / REF_LEN)
    a = beta
    y = x
    while len(y) < N:
        y = np.append(y, a*x)
        a *= beta
        
    return y[0:N+1]

y_np = KS_with_np(b, Fs*2, REF_LEN=50)
IPython.display.Audio(y_np, rate=Fs)

y_np = KS_with_np(c, Fs*2, REF_LEN=150)
IPython.display.Audio(y_np, rate=Fs)

音乐,走起！

让我们用KS算法来生成一段甲壳虫乐队著名的 opening chord of “A Hard Day’s Night”

在很多文献中，这段音乐的音为$D_3, F_3, G_3, F_4, A_4, C_5$D3​,F3​,G3​,F4​,A4​,C5​，为了让音乐听上去更“深”，我们添加了一个$D_2$D2​音

在西方音乐中，$A_4$A4​音的频率为440Hz，其他的音的频率可以根据$A_4$A4​计算得到 $f(n) = A4 \times 2^{n/12}$f(n)=A4×2n/12，其中n为目标音与A4的半音符差距个数，如果比A4高，那么n为整数，否则为负数。

我们将不同的音赋予不同的增益，就能生成一段类型的音乐啦。下面的代码需要懂一些乐理知识（我不会啊！），但是大致能够明白在干什么。freq根据音的字符串得到音的频率。

def freq(note):
    # general purpose function to convert a note  in standard notation 
    #  to corresponding frequency
    if len(note) < 2 or len(note) > 3 or \
        note[0] < 'A' or note[0] > 'G':
        return 0
    if len(note) == 3:
        if note[1] == 'b':
            acc = -1
        elif note[1] == '#':
            acc = 1
        else:
            return 0
        octave = int(note[2])
    else:
        acc = 0
        octave = int(note[1])
    SEMITONES = {'A': 0, 'B': 2, 'C': -9, 'D': -7, 'E': -5, 'F': -4, 'G': -2}
    n = 12 * (octave - 4) + SEMITONES[note[0]] + acc
    f = 440 * (2 ** (float(n) / 12.0))
    #print note, f
    return f


def ks_chord(chord, N, alpha):
    y = np.zeros(N)
    # the chord is a dictionary: pitch => gain
    for note, gain in chord.items():
        # create an initial random-filled KS buffer the note
        x = np.random.randn(int(np.round(float(Fs) / freq(note))))
        y = y + gain * KS(x, N, alpha)
    return y  

# A Hard Day's Night's chord
hdn_chord = {
    'D2' : 2.2, 
    'D3' : 3.0, 
    'F3' : 1.0, 
    'G3' : 3.2, 
    'F4' : 1.0, 
    'A4' : 1.0, 
    'C5' : 1.0, 
    'G5' : 3.5,
}
    
IPython.display.Audio(ks_chord(hdn_chord, Fs * 4, 0.995), rate=Fs)

# 另一段音乐
hdn_chord = {
    'C3' : 1.2, 
    'F#3' : 1.0, 
    'Bb3' : 1.0, 
    'E4' : 1.2, 
    'A4' : 1.0, 
    'D5' : 1.0, 
}
    
IPython.display.Audio(ks_chord(hdn_chord, Fs * 2, 0.995), rate=Fs)

重申一遍：以上所有代码和输出音频都可以在这个在线notebook进行下载和播放。