音楽プログラミングの超入門(仮)

Python / 音楽情報処理 初心者が、初心者にも分かるような記事を書きたい。

【Python】 大聖堂の響きを再現!! : インパルス応答の畳み込み

導入:インパルス応答と残響

インパルス応答(英: Impulse response)とは、インパルスと呼ばれる非常に短い信号を入力したときのシステムの出力である。インパルス反応、重み関数 (weighting function) とも。インパルスとは、時間的幅が無限小で高さが無限大のパルスである。実際のシステムではこのような信号は生成できないが、理想化としては有益な概念である。

http://ja.wikipedia.org/wiki/インパルス応答

(^^;)< ちょっと何言ってるか分かんないです。

と感じた人も少なからずいるのではないでしょうか?
ですが、インパルス応答はイメージをつかんでしまえば、それほど難しいものではありません。ここで、簡単に説明しましょう。

インパルス応答

まず、下図みたいな信号を(単位)インパルスといいます。
f:id:yukara_13:20131220092958p:plain
時間 0 のところが 1 になっているだけの信号ですね。
これをあるシステムに入力すると、下みたいに出力 h[t] が返ってくるとします。
f:id:yukara_13:20131220093033p:plain

このシステムは、時間が変わっても同じ入力に対しては同じ出力を返します(時不変)。そして、入力が2倍になったら出力も2倍に、入力が半分になったら出力も半分になります(線形)。

この二つの性質を満たすシステムを、線形時間不変なシステムといいます。
信号処理にでてくるシステムは、多分大体これだと思います。(厳密には違っても簡単に考えるためそういうことにしている場合が多いのではないでしょうか?)

そして、さっきの出力 h[t] がインパルス応答です!!
だからなんだ (^^;)・・・。と思うかもしれませんが、さらにもう少し考えてみましょう。このシステムに次のような入力が入ってきたときのことを考えます。

f:id:yukara_13:20131220093055p:plain

ここで、さっきの線形時間不変な性質を思いだしましょう。t = 1 の信号に対してはインパルス応答の2倍の出力、t = 2 の信号に対してはインパルス応答の1倍の出力が現れますね。ということは、このシステムの全体の出力は、下図のように表せますね。

f:id:yukara_13:20131220093115p:plain

(あれ?このシステムの出力、インパルス応答だけで計算できてない?)
まさにその通りですね。インパルス応答はシステムの全てを表しています(ただし線形時間不変なシステムに限る)。

まとめると、信号 x[t] をシステムに入力したとき、システムのインパルス応答を h[t] とすると、システムの出力 y[t] は次の式で計算できます。


\begin{align}
\displaystyle
 y[t] = \sum_{k=-\inf}^{\inf}x[k]h[t-k]
\end{align}

この形の演算を畳み込みといいます。
システムの出力を計算したい時は、事前にインパルス応答を測っておけばいいんですね。

残響

風呂場で歌ったりすると、とても声が響きますね。あの響きを残響といいます。
ここで風呂場を一つのシステムと考えましょう。ある音を入力すると、それに残響を付けた音を出力しますね。勘のいい人はすでに気付いているかもしれませんが、風呂場のインパルス応答を測っておけば、実際に風呂場にいかなくても、畳み込みで風呂場の音を再現できる気がしませんか?そのとおりですね。もちろん風呂場だけではなく、大聖堂でも東京ドームでも、インパルス応答さえ測っていれば、そこでの音を再現できるのです。

 では、どのようにインパルス応答を測るのでしょうか?さきほどの話にのっとると、(単位)インパルスの音を鳴らしてそれを録音すればいいですね。まあそうなんですが、理論上インパルスは長さがないので音がないです、というか生成不可能です・・。なので、近似的にインパルス応答を測ることになりますが、そのへんの詳細は省きます!!正直詳しく知らない・・ まあ、インパルスの音は「パン!」っていう何かが破裂したような音だと思ってください。

 とりあえず実際にインパルス応答だけの音を聴いてみましょう(インパルス応答は長さのある信号なので音として聴けます)。下記のように、チャペルとかのインパルス応答を配布しているサイトがあるようです。

Browse the Auralization Database | The Open Acoustic Impulse Response Library

今回は、ここから "The Lady Chapel, St Albans Cathedral" のインパルス応答をダウンロードしてみました。こんな感じの音です。

さすがチャペルといった感じで、残響がなかなか長いですね!これを適当な歌声や楽器演奏に畳み込めば、このチャペルで演奏した雰囲気が味わえるということになります。(日本にいるのに!)
というわけで、実際にやってみましょう。

音響信号の畳み込み

音響信号にインパルス応答を畳み込んでいきますが、そのまま実装すると畳み込み演算ってむちゃくちゃ重いです(式を見たら何となく分かりますよね・・)。ですが、ちゃんと高速に計算するテクニックがあるので大丈夫です(知らんけど)

 実は、時間領域での畳み込みは、周波数領域での乗算になります。なので、入力信号とインパルス応答を両方 FFT して、そいつらを複素乗算して、それを iFFT すれば割と速く畳み込みが行えます。(簡単!!)少々分かりにくかったかもしれませんが、下のサイトで超分かりやすく説明されていたので参照してください。

FIR型とインパルス応答の畳み込み : ソフトウェアと音響のデジタル信号処理

今回はこのサイトにのってる、サンプリング・リバーブというアルゴリズムで実装してみましょう。これは、音響信号全体を一気に畳み込みするのではなく、リアルタイム性とかを考慮して、フレームを少しずつずらしながら畳み込みをする方法ですね。下みたいな感じで実装してみました。

# -*- coding: utf-8 -*-
"""
Reverberation
"""
from scipy import array, ceil, float64, int16, log2, zeros
from scipy import ifft
from scipy.fftpack import fft
# from scipy.io.wavfile import read, write
from scikits.audiolab import wavread, wavwrite

# ========
#  HELPER
# ========

### Exponent of next higher power of 2
def nextpow2(n):
    return ceil(log2(abs(n)))

# =====================
#  サンプリングリバーブ
# =====================
def sampling_reverb(signal_array, impulse_response):
    sigL = len(signal_array)
    irL = len(impulse_response)

    ### インパルス応答の長さを調整
    new_irL = int((2 ** nextpow2(irL)) * 2) # 2フレーム分
    frameL = new_irL / 2
    new_IR = zeros(new_irL, dtype = float64)
    new_IR[: irL] = impulse_response

    ### 入力信号をいい感じの長さにする
    frame_num = int(ceil((sigL + frameL) / float(frameL)))
    new_sigL = frameL * frame_num
    new_sig = zeros(new_sigL, dtype = float64)
    new_sig[frameL : frameL + sigL] = signal_array

    ### インパルス応答を畳み込む
    ret = zeros(new_sigL - frame_num, dtype = float64) # 残響を付けた信号を入れていく
    ir_fft = fft(new_IR) # インパルス応答の FFT
    for ii in xrange(frame_num - 1):
        s_ind = frameL * ii
        e_ind = frameL * (ii + 2)

        sig_fft = fft(new_sig[s_ind : e_ind]) # 信号の FFT
        ### 畳み込み
        ret[s_ind : s_ind + frameL] = ifft(sig_fft * ir_fft)[frameL :].real
        

    print new_irL, sigL, new_sigL
    return ret[: sigL]


def test():
    wav_file = "../wav/harmonica.wav"
    impulse_response_file = "../data/impulseResponse/stalbans_a_short_44100.wav"

    data, fs_sig, fmt = wavread(wav_file)
    impulse_response, fs_ir, fmt2 = wavread(impulse_response_file)
    assert fs_sig == fs_ir, ":("

    ### 残響付加
    reverbed_signal = sampling_reverb(data, impulse_response)

    ### クリッピング防止
    if max(reverbed_signal) > 1:
        reverbed_signal /= (max(reverbed_signal) * 1.2)

    ### WAVファイルに書き出し
    wavwrite(reverbed_signal, "../wav/harmonica_reverb.wav", fs = fs_sig)

if __name__ == "__main__":
    test()

音を聴いてみよう

実際に上のコードで残響付加をやってみました

これが↓

適当にハーモニカを吹いたやつ

こうなる!↓

チャペルの残響を付けたやつ

どうでしょうか・・?なんだか壮大な感じになっていないでしょうか??未熟な演奏が、イキなアレンジに聴こえてこないでしょうか??小宇宙(コスモ)を感じな(ry

どう感じるかはともかく、色々な音とインパルス応答で試してみると面白いと思います。是非やってみてください。

まとめ

ただ単にエフェクトとかで残響を付けるのではなく、特定環境(好きな場所)の音響を再現できるのが面白いところですね。