オールパスフィルターによるリバーブの実装

本記事では、オールパスフィルタによるリバーブを実装してみました。オールパスフィルタとは振幅特性を変えずに位相特性を変えるフィルタのことです。

くし型フィルタによるリバーブの実装を以前紹介しましたが、シュレーダーの論文 [1] ではくし型フィルタの欠点をいくつか指摘しており、その欠点を改善する方法としてオールパスフィルタによるリバーブを提案しています。

この記事では、文献 [1] で記載されていた「くし型フィルタの欠点」と「オールパスフィルタによるリバーブ」について紹介していきます。

1 くし型フィルタの欠点
- 1.1 リバーブの振幅特性が平坦ではない
- 1.2 エコー密度が小さい
2 オールパスフィルタ
- 2.1 オールパスフィルタについて
- 2.2 エコー密度を高める方法
3 プログラム
- 3.1 ソースコード
- 3.2 実行方法
4 処理結果
5 おわりに

くし型フィルタの欠点

文献 [1] ではくし型フィルタの欠点を以下の２つ指摘しています。

フィルタの振幅特性が平坦ではない
エコー密度が小さい

エコー密度とは、1つのパルスが入力されたときのリバーブエフェクタが出力する1秒あたりの反射音の数です。

リバーブの振幅特性が平坦ではない

くし型フィルタである以下の式

$$
y[n] = x[n] + g y[n-D] $$

は下図のようにくし型の振幅特性をしています。

上図のようにくし型フィルタの振幅特性は平坦ではないため、カラレーションのように音色が変化してしまう欠点があります。

この音色の変化をなくすために、振幅特性が平坦なフィルタを使う必要があります。

補足：私も知りませんでしたが、論文中に出てきたカラレーション (coloration) とは以下のような現象みたいです。

カラレーション

室内で起こる音響現象の一つです。

発せられた直接音と壁などによる反射音が干渉して、音色が変化して聞こえる現象です。中・高音域で多重反射が起きている室内などで発生しやすい。

引用元：【カラレーション】 - オーディオ辞典

エコー密度が小さい

くし型フィルタのインパルス応答を以下に示す。

上図のようにサンプル番号Ｄごとに値があるようなインパルス応答となっています。

遅延時間$\ \tau$ を 40ms （サンプリング周波数 48kHzの場合、D=48x40=1920）に設定した場合、エコー密度は1000ms/40ms=25です。このエコー密度は実際の残響に比べて小さく、音がはためいているように聞こえてしまいます。

文献[1] によれば、はためきをなくすにはエコー密度が約1000必要であり、1000あればそれ以上のエコー密度のリバーブと聴き分けができないみたいです。そのため、リバーブの改善にはエコー密度を1000に近づける必要があります。

補足：「音がはためいている」は文献 [1] の"fluttering"の直訳です。

オールパスフィルタ

オールパスフィルタとオールパスフィルタでエコー密度を高める方法を説明したいと思います。

オールパスフィルタについて

オールパスフィルタとは冒頭で説明したように周波数の振幅特性を変えずに位相特性を変えるフィルタのことです。

くし型フィルタの式の直接音 x[n] を$-g $ 減衰させ、 D だけ遅延させた直接音 x[n-D] を加えることでオールパスフィルタは作成できます。

$$
\begin{align}
●& くし型フィルタ \\[5pt] &y[n] = x[n] + g y[n-D] \\[9pt] ●& オールパスフィルタ \\[5pt] &y[n] = \color{red}{-g x[n]} + \color{red}{x[n-D]} + g y[n-D] \\
\end{align}
$$

オールパスフィルタのブロック図とインパルス応答は以下です。

z 変換した式は以下です。

$$
H(z) = \frac{-g+z^{-D}}{1-g z^{-D}}
$$

$z=e^{jω}$ を代入して、絶対値を計算すると以下のように振幅特性が得られます。

$$
\begin{align}
|H(\omega)| &= \frac{|-g+e^{-j\omega D}|}{|1-g e^{-j\omega D}|} \\[5pt] &= |e^{-j\omega D}|\frac{|1-g e^{j\omega D}|}{|1-g e^{-j\omega D}|} \\[5pt] &= 1 \cdot \frac{\sqrt{(1-g \cos{(\omega D)})^2+(g \sin{(\omega D)})^2}}{\sqrt{(1-g \cos{(\omega D)})^2+(g \sin{(\omega D)})^2}} \\[5pt] &= 1
\end{align}
$$

上記のように振幅特性がすべての周波数で1となっているため、このオールパスフィルタをリバーブ作成に使用します。

エコー密度を高める方法

エコー密度を高めるためにオールパスフィルタを以下のように縦続接続します。

各オールパスフィルタの遅延時間 D2, D3, …については前の段の遅延時間の約1/3に設定します。つまり、n 段目のオールパスフィルタの遅延時間 Dn は数式で以下のように表せます。

$$
D_{n} = D_{1}\left(\frac{1}{3}\right)^{n-1} \tag{A}
$$

ただし、反射音の打ち消しや重なりが起こらないように各段の遅延時間は厳密に(A)式に従わず、微妙にずらして設定します。

このとき、エコー密度 d は以下のように表せます。

$$
d \approx \frac{3^{n-1}}{\tau_1}
$$

$\tau_1=100$ ms、オールパスフィルタを5段接続した場合、エコー密度は約810となり、1000に十分近づきます。

各減衰率 $g_n$ については 0.7 が良いと記載がありました。

補足：各段の遅延時間どうしは反射音の打ち消しが起こらないようにそれらの最小公倍数をなるべく大きくしたほうがよいです。一番ベストなのは互いに素となる遅延時間を使うことです。

例として、D1=4800、g=0.7、オールパスフィルタを２つ接続したときのインパルス応答を以下に示します。上図がD2=D1/3のとき、下図がD2=D1/3+1のときのインパルス応答です。

図：遅延時間を微妙にずらす理由 — 上図：D2 = D1/3 の例　下図：D2 = D1/3+1 の例

図からわかるように微妙にずらしていない場合、反射音が打ち消しあい、エコー密度が小さくなっています。特に赤く示したところは顕著に小さくなっています。

このように、反射音が打ち消しあい、エコー密度も小さくなってしまいますので、遅延時間は(A)式から微妙にずらします。

プログラム

オールパスフィルタによるリバーブを適用するプログラムをPythonで実装しました。ソースコードと実行方法について説明します。

ソースコード

入力データにオールパスフィルタを5つ縦続接続したリバーブを適用するプログラムは以下です。

import soundfile as sf
import numpy as np
import scipy.signal as sg

# パラメータ
tau1 = 100    # 1段目遅延時間 [ms]
gain = 0.7    # 減衰率
unit_num = 5  # 縦続接続数
wav_in_name  = "input_short.wav"  # 入力音データ名
wav_out_name = "output_short.wav" # 出力音データ名

# WAVファイルを読み込む
x, fs = sf.read(wav_in_name)

# [ms] → [サンプル] に変換
D = [0] * unit_num
D[0] = int(tau1 * fs / 1000)
for i in range(1, unit_num):
    D[i] = D[i-1]//3+1  # 重ならないように+1してずらす

# フィルタ係数を設定
for i in range(unit_num):
    b = np.zeros(D[i]+1)
    a = np.zeros(D[i]+1)
    b[0] = -1.0*gain
    a[0] = 1.0
    a[D[i]] = -1.0*gain
    b[D[i]] = 1.0
    x = sg.lfilter(b, a, x) # リバーブをかける

# wavファイルを書き込む
y = x    
y = y/np.max(np.abs(y))   # ノーマライズ
sf.write(wav_out_name, y, fs, subtype='PCM_16')

18～19行目：2段目、3段目...の遅延時間 Dn を設定します。前の段の遅延時間の1/3に+1して、反射音の重なりが起こらないようにします。

21～29行目：オールパスフィルタの接続数だけフィルタをかけています。

33行目：オーバーフローしたときのためにノーマライズしています。

実行方法

(1) プログラムを実行するディレクトリにソースコード(allpass_reverb.py)と入力 WAV データを格納する。

(2) ソースコード5～10行目のパラメータと入出力データ名を修正する。

# パラメータ
tau1 = 100    # 1段目遅延時間 [ms]
gain = 0.7    # 減衰率
unit_num = 5  # 縦続接続数
wav_in_name  = "input_short.wav"  # 入力音データ名
wav_out_name = "output_short.wav" # 出力音データ名

(3) 以下のコマンドで python を実行することで、リバーブをかけたWAVデータが出力される。