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

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

什么是Karplus-Strong算法

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

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

import numpy as np
import matplotlib.pyplot as plt
import IPython
%matplotlib inlinedef KS_1(x, N):'''x: 初始序列N: 输出序列的最大长度'''y = xwhile len(y) < N:y = np.append(y, x)y = y[0:N+1]return yFs = 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] y[n]=x[n]+y[n−M]
并假设x是 finite-support 信号
x[n]={0n<0bn0≤n<M0n≥Mx[n] = \begin{cases} 0 & n< 0\\ b_n & 0 \le n < M \\ 0 & n \ge M\\ \end{cases} x[n]=⎩⎪⎨⎪⎧0bn0n<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]=α⌊n/M⌋x[nmod  M]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，求得betay = 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 = betay = xwhile len(y) < N:y = np.append(y, a*x)a *= betareturn 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”

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

在西方音乐中，A4A_4A4音的频率为440Hz，其他的音的频率可以根据A4A_4A4计算得到 f(n)=A4×2n/12f(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 frequencyif len(note) < 2 or len(note) > 3 or \note[0] < 'A' or note[0] > 'G':return 0if len(note) == 3:if note[1] == 'b':acc = -1elif note[1] == '#':acc = 1else:return 0octave = int(note[2])else:acc = 0octave = 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]] + accf = 440 * (2 ** (float(n) / 12.0))#print note, freturn fdef ks_chord(chord, N, alpha):y = np.zeros(N)# the chord is a dictionary: pitch => gainfor note, gain in chord.items():# create an initial random-filled KS buffer the notex = 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进行下载和播放。