６．コンサートホールの反響を再現しよう - keep learning blog（キープラーニングブログ）

――音楽だけが世界語であり、翻訳される必要がない。そこにおいては魂が魂に話しかける。――

コンサートホールの反響を再現しよう

コンサートホールの反響を再現しよう

第４回の記事で紹介したとおり、私はいま音響学の教科書を一冊読み終えて、勝手にスクリプト（C言語等と違ってコンパイル不要なプログラムのこと）を作って遊んで自習しています。

プログラミングは素人同然なので、しょせんは趣味レベルではありますが、このブログでもときどき家で書いたソースコードを公開していきたいと思っています。今回はその第一弾です。

f:id:yuki0718:20190728045010j:plain

環境構築

私がプライベートで使用しているノートPCは、社会人3年目くらいに買った、東芝 dynabookのKIRA V63/27Mという古めの機種です。容量の大きいSSDに換装したり、Windows10を買ってOSを入れ直したりしながら、長く愛用しています。

一応CPUはCORE i5が入っているので、Pythonでちょっとしたスクリプトを動かしてもぎりぎり堪えます。私はMacを使ったことがないのですが、「Windowsマシンだけど見た目がMacbookっぽくてオシャレ」と人気の機種でした（過去形）。

kakaku.com

さて、私のような軟派なWindowsユーザーでも、Pythonというプログラミング言語は簡単に始めることができます。私は2年前にAnacondaという（これまた蛇っぽい）名前のディストリビューションを使って、Pythonを動かすための環境を構築しました。

環境構築と聞いてピンとこない方もいらっしゃると思いますが、プログラミング言語というのは、すべての処理を一から自分で書くものではありません。むしろ、開発者さんたち優れた先人の作った環境を間借りして、その環境の中で自分が興味のある処理を命令し、実行します。

これはちょうど、左官屋さんが大工さん（開発者）の作った骨組みにモルタルで壁を塗るような作業です。左官屋さんは一から家を作りません。しかし、大工さんが骨組みを作ってくれているので、好みのデザインで表面を塗装して、堂々と自分の仕事を人に見せることができます*1。

環境構築の方法については、私などが説明するまでもなく、説明の上手なプロの方々が非ITエンジニア向けに分かりやすく書いてくださっているので、ぜひ参考にしてみてください。

tonari-it.com

programming-study.com

今回は、私の読んだ「Introduction to Audio Processing」の中で紹介されていたアルゴリズムのうち、Pythonで簡単に実装できそうな具体例として、Schroeder's reverb（シュレーダー・リバーブ）をご紹介します。

Schroeder's reverb

反響（Reverberation）とは、コンサートホールのような広い閉空間で演奏を行ったときに発生する空間音響効果の一つです。天井や壁で何度も跳ね返った音の波が、観客の耳に届いたときのなんともいえない響きのことです。

歴史的にも、反響を人工的に生み出そうとする研究は数多く存在しています。その中でも特に有名なのが、1962年にManfred R. Schroederさんの発表した論文「Natural Sounding Artificial Reverberation」の中で示された、Schroeder's Reverbと呼ばれるアルゴリズムです。

Schroeder's Reverbは、4つの並列なくし形フィルタ（Comb Filter）と2つの全透過フィルタ（All-Pass Filter）とを直列に組み合わせたデジタルフィルタを通過させるだけで、まるでコンサートホールにいるような反響が付加されるというアイディアです。

一般に、デジタルフィルタが単純な多項式（Polynomial）の差分方程式で表現される場合、出力 $y_{n}$ は、入力 $x_{n}$ を使って次のように書けます。

一般的なフィルタの差分方程式 $y_{n} - a_{1}y_{n-1} - a_{2}y_{n-2} \cdots -a_{i}y_{n-i} = x_{n} + b_{1}x_{n-1} + b_{2}x_{n-2} \cdots +b_{j}x_{n-i}$

$n$ は時間を表す整数です。コンピュータは連続量を入力できないので、時間軸は離散的（デジタル）になります。 $a_{i}$ や $b_{j}$ というのはフィルタの特性を決定する係数です。 $y_{n}=Y(ω)e^{iωn}$ と $x_{n}=X(ω)e^{iωn}$ を代入して周波数（音の高さω）の空間に移行させると、出力と入力の比で定義されるフィルタ応答関数 $H(ω) = Y(ω)/X(ω)$ の分母は $a_{i}$ のみで書け、分子は $b_{j}$ のみで書けることが分かります。

さて、上述したくし形フィルタというのは、ド⇒ミ⇒ソ⇒シ⇒・・・といった一定の間隔で、櫛のようなピークが並んだフィルタのことです。上述した分母が0に近づく（発散する）ところで櫛ピークが出るので、 $a_{i}$ の値（以下に示すG,Dの値）を上手に設計することが肝要です。

くし形フィルタの差分方程式 $y_{n} - Gy_{n-D} = x_{n}　(G, Dは定数)$

f:id:yuki0718:20190804040936p:plain — くし形フィルタの応答関数（周波数特性）

他方、全透過フィルタというのは、すべての高さの音について透過するものの、高さの異なる波同士の位相を微妙にずらす作用をするフィルタです。音のハーモニーのバランスをうまい具合に崩す、というイメージです。

全透過フィルタの差分方程式 $y_{n} - Gy_{n-D} = -Gx_{n} + x_{n-D}　(G, Dは定数)$

フィルタの設計は非常に難しい問題ですが、このG、Dといった定数の値をうまい具合に決定する方法が、Schroederさんの論文に書かれているので、そのとおりに実装すれば誰でも簡単に反響効果を得ることができます。

PythonでSchroeder's reverbの実装

いきなり結論ですが、ソースコード載せます。すごく簡単なので、Pythonのコードについて特に説明はしません（なんだそれ）。私はAnacondaに入っている「jupyter notebook」というエディタで書いています。書いて実行、書いて実行、というテストに適したエディタですので、大変おすすめです。

なお、このスクリプトはsoundfileという音楽読込み用ライブラリと、numpyという数値計算用ライブラリと、scipyのsignalという信号処理用ライブラリをAnacondaでインストールしないと動きません*2。

import soundfile as sf
import IPython.display as ipd
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal as sg

#■音響データにシュレーダー・リバーブを付加する関数を定義■
def Schroeder_Reverb(wavdata, Fs):
    
    #Combフィルターを4つ定義する
    D1 = 1757 #39.85msec
    G1 = -0.871402
    C1 = 0.25
    #フィルタの応答関数(H(z)=Σbk*z^-k/Σak*z^-k)を分子b分母aに分けて設定
    b1 = [1]
    a1 = [1]
    a1.extend([0]*(D1-1))
    a1.extend([G1])
    
    #残りの3つも同様にパラメータを設定
    D2 = 1592 #36.10msec
    G2 = -0.882762
    C2 = 0.25
    b2 = [1]
    a2 = [1]
    a2.extend([0]*(D2-1))
    a2.extend([G2])
    
    D3 = 1467 #33.27msec
    G3 = -0.891443
    C3 = 0.25
    b3 = [1]
    a3 = [1]
    a3.extend([0]*(D3-1))
    a3.extend([G3])
    
    D4 = 1330 #30.15msec
    G4 = -0.901117
    C4 = 0.25
    b4 = [1]
    a4 = [1]
    a4.extend([0]*(D4-1))
    a4.extend([G4])
    
    #並列にCombフィルタをかける(応答関数の分子と分母をリスト形式でlfilter関数に渡す)
    output1 = C1 * sg.lfilter(b1, a1, wavdata)
    output2 = C2 * sg.lfilter(b2, a2, wavdata)
    output3 = C3 * sg.lfilter(b3, a3, wavdata)
    output4 = C4 * sg.lfilter(b4, a4, wavdata)
    output = output1 + output2 + output3 + output4
    
    #続けて2つのAll-Pass Filterのパラメータも定義する
    D5 = 220; #5.0msec
    G5 = -0.7
    #フィルタの応答関数(H(z)=Σbk*z^-k/Σak*z^-k)を分子b分母aに分けて設定
    b5 = [G5]
    b5.extend([0]*(D5-1))
    b5.extend([1])
    a5 = [1]
    a5.extend([0]*(D5-1))
    a5.extend([G5])
    
    D6 = 75 #1.7msec
    G6 = -0.7
    b6 = [G6]
    b6.extend([0]*(D6-1))
    b6.extend([1])
    a6 = [1]
    a6.extend([0]*(D6-1))
    a6.extend([G6])
    
    #直列にAll-Passフィルタをかけて戻り値に出力
    output = sg.lfilter(b5, a5, output)
    return sg.lfilter(b6, a6, output)

#■メイン処理■
if __name__ == "__main__":
    
    #ファイルパス
    fname = './guitar.wav'
    
    #波形データwavdataとサンプリング周波数Fsを取得
    data, Fs = sf.read(fname)
    #ステレオ音源の場合はチャンネル1を取得
    if(data.shape[1] == 1):
        wavdata = data
    else:
        wavdata = data[:, 0]
    
    #IPythonライブラリを使ってjupyter notebook上に再生バーを表示
    print('Original Sound')
    #ipd.display(ipd.Audio(data=wavdata, rate=Fs))
    
    #横軸用の時間を作っておく
    time = np.arange(0, len(wavdata))/Fs
    
    #matplotというライブラリを使って音源の波形グラフを表示
    plt.figure(figsize=(16, 5))
    plt.plot(time, wavdata)
    plt.title('Wave data of the original music', fontsize=12)
    plt.xlabel('Time [sec]')
    plt.ylabel('Amplitude')
    plt.xlim(0, 5)
    plt.show()
    
    #短時間フーリエ変換(STFT)を実行
    #窓関数の種類はハニング窓,窓の大きさは20msec,窓の重なり幅は10mecで設定
    F, T, Adft = sg.stft(wavdata, fs=Fs, window='hann', nperseg=0.02*Fs, noverlap=0.01*Fs)
    
    #振幅をdBに変換
    P = 10 * np.log(np.abs(Adft))
    
    #STFTのスペクトログラムを表示
    plt.figure(figsize=(16, 5))
    plt.title('Spectrogram of the original music', fontsize=12)
    plt.xlabel('Time [sec]')
    plt.ylabel('Frequency [kHz]')
    plt.xlim(0, 5)
    plt.pcolormesh(T, F/1000, P, cmap = 'rainbow')
    plt.colorbar(orientation='horizontal').set_label('Amplitude [dB]')
    plt.show()
    
    #自己定義関数でReverbを付加したサウンドを生成する
    reverb_wavdata = Schroeder_Reverb(wavdata, Fs)
    
    #IPythonライブラリを使ってリバーブ付加後の波形グラフを表示
    print('Reverbed sound')
    #ipd.display(ipd.Audio(data=reverb_wavdata, rate=Fs))
    
    #matplotというライブラリを使って音源の波形グラフを表示
    plt.figure(figsize=(16, 5))
    plt.plot(time, reverb_wavdata)
    plt.title('Wave data of the reverbed music', fontsize=12)
    plt.xlabel('Time [sec]')
    plt.ylabel('Amplitude')
    plt.xlim(0, 5)
    plt.show()
    
    #短時間フーリエ変換(STFT)を実行
    #窓関数の種類はハニング窓,窓の大きさは20msec,窓の重なり幅は10mecで設定
    F2, T2, Adft2 = sg.stft(reverb_wavdata, fs=Fs, window='hann', nperseg=0.02*Fs, noverlap=0.01*Fs)
    
    #振幅をdBに変換
    P2 = 10 * np.log(np.abs(Adft2))
    
    #STFTのスペクトログラムを表示
    plt.figure(figsize=(16, 5))
    plt.title('Spectrogram of the reverbed music', fontsize=12)
    plt.xlabel('Time [sec]')
    plt.ylabel('Frequency [kHz]')
    plt.xlim(0, 5)
    plt.pcolormesh(T2, F2/1000, P2, cmap = 'rainbow')
    plt.colorbar(orientation='horizontal').set_label('Amplitude [dB]')
    plt.show()

実行結果

実行結果はこちらです。まずはオリジナルの音声（フリー素材）がこちら。

そしてSchroeder's reverbのフィルタを通過して出てきた音声がこちら。

うまく言葉にできませんが、なかなかのコンサートホール感が出ていると思います。先ほどのフィルタの係数をうまく調整することで、様々な広さのホールにいるような雰囲気を生み出すことができます。

また、このスクリプトでは、フィルタを通る前後で音の波形がどう変化するかも出力されます。まずはオリジナルの時間‐振幅グラフと、波を0.02秒毎に区切ってフーリエ変換（短時間フーリエ変換STFT）したものを並べた、いわゆるスぺクトログラムのグラフを示します。

f:id:yuki0718:20190728035921p:plain

f:id:yuki0718:20190728035947p:plain

ギターを細かにかき鳴らしている演奏なので*3、弦が弾かれたところが赤くなっています。また、次回の記事で紹介しますが、ギターの音は固有周波数（一番低い音）の整数倍の周波数にピークが発生するので、縦軸方向に等間隔な横線が並んでいるのが見て取れます。

続いて、Schroeder's reverbのフィルタを通過して出てきた音声のグラフを示します。

f:id:yuki0718:20190728040738p:plain

f:id:yuki0718:20190728040801p:plain

波形からして、明らかに音に厚みが出ているのが分かると思います。コンサートホールの天井や壁で反射したときのように、音の波が重なり合っています。スペクトログラムも横の線が多くなって、かなり重層的な音声データという印象になりました。

以上、誰でもかんたんに実装できる、Pythonで音声処理その１でした。

大学の研究は「Matlab」という学術用の有料ソフトを使ってやっているので、家では使えません。音声処理の雰囲気を感じてもらえればと思い、家でも使えるPythonで書いてみました。また、ブログで数式やコードをどこまで載せられるか、テストを兼ねています（きれいに掲載できてますかね？）。

ちょっと専門的な話になってついて行けなかった方はすみません。プログラミングの筆がのってきたので、次回はシンセサイザーのようにギター風の音色を生成するアルゴリズムを紹介したいと思います。

稚文をお読みいただきありがとうございました。

*1:自分の親戚に大工や建築家がいるので喩えに使ってみました。私もそのくらい器用でクリエイティブな人間なら良かったのですが…。

*2:IPythonのdisplayというライブラリは、単に私が使っているjupyter notebookというエディタ上に音声の再生ボタンを作るためなので動作には関係ありません。

*3:私の弟は高校から軽音部でエレキギターを弾いており、実家に帰るとものすごい超絶技巧を高速で奏でています。ユーロ圏のヘヴィメタルが好きみたいですが、私には良さがいまいち分かりません。