幅と中心を個別に調整できるバンドリジェクトフィルターの作り方

中心コントロールは、減衰が最大になる周波数を設定します。幅コントロールは、慣例的に−3 dB点と呼ばれる2つの半電力境界の間隔を設定します。ノッチはバンドストップフィルターの狭帯域版です。SciPyの公式ドキュメントでは、品質係数は中心周波数をこの帯域幅で割った値と定義されています。SciPy: iirnotch

コントロールを設計する前に、幅の単位を選びます。ヘルツで表す幅は絶対的な周波数範囲を保ちます。オクターブで表す幅は周波数比を示します。W3CのAudio EQ Cookbookは、ノッチフィルターのオクターブ帯域幅に対応しており、デジタル係数の式で周波数ワーピングを考慮しています。これらは動作の異なるコントロールなので、選んだ単位を明確に表示してください。W3C: Audio EQ Cookbook

ここではヘルツを使います。中心周波数が1,000 Hz、幅が200 Hzなら、対応するQは5です。幅を200 Hzに保ったまま中心周波数を2,000 Hzに移動するには、Qを10にする必要があります。Qを5に固定すると、帯域幅は代わりに400 Hzになります。この例の計算から、Q固定のコントロールでは今回選んだ要件を満たせないことがわかります。

2次の無限インパルス応答フィルター、別名バイクアッドを使います。分子係数が3つ、分母係数が3つあります。最初の分母係数を1に正規化すると、サンプルごとの式は次のようになります。

[ y[n]=b_0x[n]+b_1x[n-1]+b_2x[n-2]-a_1y[n-1]-a_2y[n-2]. ]

これはAudio EQ Cookbookに記載されているバイクアッドの関係式です。W3C: Audio EQ Cookbook

サンプルレート \(F_s\)、中心周波数 \(f_c\)、幅 \(B\) に対して、次を計算します。

\[ \omega_c=2\pi f_c/F_s,\qquad \beta=\tan(\pi B/F_s),\qquad g=\frac{1}{1+\beta}. \]

次に、次の値を使います。

\[ b=g[1,-2\cos(\omega_c),1], \]

\[ a=[1,-2g\cos(\omega_c),2g-1]. \]

これらはSciPyの実装で使われているノッチ係数の式を、帯域幅を直接入力する形に書き換えたものです。SciPy filter-design source

どちらのコントロールを変えても、複数の係数が変わる場合があります。独立性とは、結果として得られる応答に関するものです。幅を変えても減衰周波数は一定に保たれ、中心周波数を変えても指定した帯域幅が保たれます。各ノブを1つの係数だけに割り当てても、この動作は実現できません。

次の完全なプログラムを band_reject.py として保存します。Python 3.8以降が必要で、追加パッケージは不要です。サンプルレートは48,000 Hzで、生成される入力には300、1,000、5,000 Hzの3つのトーンが含まれます。各トーンの振幅は0.2です。

このプログラムはフィルター応答を直接評価し、二分法で2つの半電力境界を求め、生成した音声を2秒間処理します。WAVの書き出しには、Pythonのドキュメントに記載された wave インターフェースを使い、モノラルの16ビットPCMファイルを作成します。出力ファイルがすでに存在する場合は、上書きせずにエラーになります。Python: wave — WAVファイルの読み書き

import argparse
import cmath
import math
import struct
import wave


def coefficients(center, width, fs):
    if not all(math.isfinite(v) for v in (center, width, fs)):
        raise ValueError('Parameters must be finite')
    if not (fs > 0 and 0 < center < fs / 2 and 0 < width < fs / 2):
        raise ValueError('Require fs > 0 and center, width inside (0, fs/2)')
    beta = math.tan(math.pi * width / fs)
    g = 1 / (1 + beta)
    c = math.cos(2 * math.pi * center / fs)
    return (g, -2 * g * c, g), (1.0, -2 * g * c, 2 * g - 1)


def magnitude(f, b, a, fs):
    z = cmath.exp(-2j * math.pi * f / fs)
    return abs((b[0] + b[1]*z + b[2]*z*z) /
               (a[0] + a[1]*z + a[2]*z*z))


def edge(lo, hi, b, a, fs):
    target = 1 / math.sqrt(2)
    for _ in range(60):
        mid = (lo + hi) / 2
        if ((magnitude(lo, b, a, fs) - target) *
                (magnitude(mid, b, a, fs) - target)) > 0:
            lo = mid
        else:
            hi = mid
    return (lo + hi) / 2


def process(samples, b, a):
    x1 = x2 = y1 = y2 = 0.0
    result = []
    for x in samples:
        y = b[0]*x + b[1]*x1 + b[2]*x2 - a[1]*y1 - a[2]*y2
        result.append(y)
        x2, x1, y2, y1 = x1, x, y1, y
    return result


def write_wav(path, samples, fs):
    if any(not math.isfinite(x) or abs(x) > 1 for x in samples):
        raise ValueError('Output exceeds PCM range; lower the input level')
    data = b''.join(struct.pack('<h', round(x * 32767)) for x in samples)
    with open(path, 'xb') as stream:
        with wave.open(stream, 'wb') as wav:
            wav.setnchannels(1)
            wav.setsampwidth(2)
            wav.setframerate(fs)
            wav.writeframes(data)


def main():
    parser = argparse.ArgumentParser()
    parser.add_argument('--center', type=float, default=1000)
    parser.add_argument('--width', type=float, default=200)
    parser.add_argument('--write', action='store_true')
    args = parser.parse_args()
    fs = 48000
    b, a = coefficients(args.center, args.width, fs)
    low = edge(0, args.center, b, a, fs)
    high = edge(args.center, fs / 2, b, a, fs)
    print(f'Edges: {low:.3f}, {high:.3f} Hz; width: {high-low:.3f} Hz')
    assert math.isclose(high-low, args.width, abs_tol=1e-5)
    assert magnitude(args.center, b, a, fs) < 1e-7
    tones = (300, 1000, 5000)
    for f in tones:
        db = 20 * math.log10(max(magnitude(f, b, a, fs), 1e-12))
        print(f'{f} Hz: {db:.2f} dB')
    samples = [sum(0.2 * math.sin(2*math.pi*f*n/fs) for f in tones)
               for n in range(2*fs)]
    filtered = process(samples, b, a)
    assert len(filtered) == len(samples)
    assert all(math.isfinite(y) for y in filtered)
    if args.write:
        write_wav('input-demo.wav', samples, fs)
        write_wav('filtered-demo.wav', filtered, fs)


if __name__ == '__main__':
    main()

応答を確認するには、python band_reject.py --center 1000 --width 200 を実行します。--write を追加すると、input-demo.wav と filtered-demo.wav が作成されます。これらのファイルを別々のトラックとして音声エディターに読み込み、開始位置を揃えて、再生ゲインを変えずに比較します。

ファイルを書き出さずにプログラム全体を実行しました。帯域幅、中心周波数での減衰、出力長、有限値出力のアサーションはすべて通りました。この実行ではWAV書き出しは試していません。表示された応答は、300 Hzで約−0.02 dB、5,000 Hzで約−0.01 dBでした。中心値は表示下限の−240 dBに達しました。この下限は数値表示上の限界であり、音声ファイルの減衰深度を測定した値ではありません。

基準設定、幅を変えずに中心周波数を変えた設定、中心周波数を変えずに幅を変えた設定の3つをテストします。これにより各コントロールの動作を個別に確認できます。次の値は、48 kHzで係数計算ルーチンと境界検索ルーチンを実行して算出したものです。録音の測定値ではなく、数学的な応答を示しています。

中心周波数の設定	幅の設定	下側の半電力境界	上側の半電力境界	算出された幅
1,000 Hz	200 Hz	904.959 Hz	1,104.959 Hz	200.000 Hz
2,000 Hz	200 Hz	1,902.441 Hz	2,102.441 Hz	200.000 Hz
1,000 Hz	400 Hz	819.692 Hz	1,219.692 Hz	400.000 Hz

境界が正確には center ± width/2 にならない点に注目してください。ここでは「中心」は減衰のヌル点を示し、「幅」は半電力交点の間隔を示します。算出された応答から、その算術的な中点がヌル周波数と一致するとは限らないことがわかります。幅の半分を単純に加減して、境界を表示しないでください。

SciPyが使える環境では、返された係数を scipy.signal.freqz で評価すると、さらに数値的な確認ができます。fs と明示的な周波数を指定し、正確な中心周波数も含めてください。ドキュメントによれば、明示的な周波数配列を指定すると、その周波数で応答が評価されます。粗い既定のグリッドでは、狭いヌル点を見逃すことがあります。SciPy: freqz

編集対象の音から、短く再現可能な部分を使います。まず減らしたい周波数領域を特定し、中心周波数をそこに設定します。幅を狭くして始め、狙った領域がカバーされるまで徐々に広げます。変更のたびにバイパスと比較し、比較が有効になるよう再生ゲインを一定に保ちます。

次の調整を選ぶ際は、この判断表を使ってください。

比較してわかったこと	次の調整
不要なトーンが残る一方、周辺の音も変わる	幅を広げる前に中心周波数を再確認する
中心周波数は適切だが、不要な領域がその外側まで広がっている	幅を少しずつ広げる
周囲の音色が消えすぎる	幅を狭めて再比較する
離れた複数のトーンを減らす必要がある	それぞれの周波数に個別のノッチ段を使うことを検討する
広く平坦な阻止領域が必要である	フィルター形状と次数を見直す

これらは実際の録音から得られた測定値ではなく、実用的な調整の提案です。3トーンの例は意図的に単純にしてあります。複雑な音について判断する前に、プロセッサが正しく動作することを確認できます。

単一のノッチには、個別に調整可能な深さコントロールもありません。Cookbookには、ノッチフィルターとピーキングEQフィルターそれぞれの係数群が示されています。ヌル点を作るのではなく、調整可能で適度なカットが必要な場合は、ピーキングEQも検討に適した設計です。W3C: Audio EQ Cookbook

同じパラメーター仕様を維持します。center_hz と width_hz を一緒に係数関数へ渡します。どちらかのパラメーターが変わるたびに係数一式を更新し、処理時のサンプルレートが変わった場合も再計算します。このプロトタイプでは、どちらのコントロールも0より大きく、サンプルレートの半分より小さい値にする必要があります。無効な設定では例外が発生します。

付属の process() 関数は、モノラルのバッファ全体を処理し、呼び出しごとに履歴をゼロから開始します。ストリーミング版では、ブロック間で遅延履歴を保持する必要があります。SciPyの lfilter には、この用途の初期遅延状態と返される最終状態が記載されています。音声チャンネルごとに別々の履歴を管理してください。SciPy: lfilter

コントロールを動かす処理は、追加の実装工程として扱ってください。この例で検証するのは固定設定であり、連続的なオートメーションは実装していません。実用的な切り替え方法の1つは、旧フィルターと新しい固定フィルターを短時間並行して動かし、新しいフィルターで履歴が蓄積されてから出力をクロスフェードすることです。実際の素材と想定する操作速度で切り替えを評価してください。

複数のフィルターを直列に接続する場合は、カスケードを1つの高次多項式に展開せず、それぞれを独立した2次セクションとして保ちます。SciPyは、数値上の問題が少ないことから、ほとんどのフィルタリング用途で2次セクションを推奨しています。SciPy: lfilter

ノッチが誤った周波数に現れる場合は、係数計算に使ったサンプルレートと音声の実際のサンプルレートが一致しているか確認します。中心周波数を動かすと幅が広がる場合は、ラッパーがQを固定していないか確認します。応答表示が浅く見える場合は、設計を変更する前に正確な中心周波数で評価します。

バッファを繰り返し処理すると不連続が発生する場合は、状態の扱いを確認してください。このデモでは process() を呼び出すたびに履歴をリセットします。WAV書き出しで範囲エラーが発生した場合は、入力振幅を下げて再実行します。このプログラムは、範囲外のサンプルを黙ってクリッピングすることを意図的に拒否します。別のファイルを書き出すには、既存のデモファイルの名前を変更するか削除する必要があります。

中心周波数のテストで指定した位置にヌル点が現れ、2つの境界交点の距離が指定した幅に保たれ、処理後の音が意図した編集を実現していれば、実装は完了です。このオフラインプロトタイプから対話型の音声ツールに移行する際も、2つのコントロールとともにこれらの確認項目を維持してください。

• Audio EQ Cookbook — W3C
• iirnotch — SciPyマニュアル
• _filter_design.py — SciPy v1.18.0ソース
• freqz — SciPyマニュアル
• lfilter — SciPyマニュアル
• wave — WAVファイルの読み書き — Pythonドキュメント