同じIIRフィルタなのに、ブロック境界で波形が崩れる理由
センサーや音声のデータを128サンプルずつ受け取り、各ブロックに signal.sosfilt(sos, block) を適用すると、境界で波形が不自然に変わることがあります。係数の設計が正しくても、過去の入力が残したフィルタ状態を毎回捨てていると連続した処理になりません。
本記事はフィルタの設計式よりも、設計済みIIRフィルタを途切れず動かす方法に絞ります。基本的な次数・遮断周波数の設計は バターワースフィルタ 、通過帯域の選択は バンドパスフィルタ を参照してください。
SOSの係数と状態は別のもの
SOS（Second-Order Sections）は、高次のフィルタを2次セクションの直列接続で表したものです。各行には分子と分母の係数が並びます。
\[ H(z) = \prod_{k=1}^{K} \frac{b_{0,k}+b_{1,k}z^{-1}+b_{2,k}z^{-2}} {1+a_{1,k}z^{-1}+a_{2,k}z^{-2}} \]ここではSciPyで設計した、分母の先頭係数が1のSOSを使います。6次ローパスなら3セクションで、係数配列は (3, 6) です。
from scipy import signal
sos = signal.butter(6, 40, fs=1000, output="sos")
print(sos.shape) # (3, 6)
係数は周波数応答を決めます。一方、状態は「直前まで何を入力したか」を記憶します。SciPyの転置直接形IIでは各セクションに2個の状態があり、1次元入力なら状態配列は (K, 2) です。SOS形式の係数を選ぶ理由は、高次の伝達関数を1つの多項式で扱うより数値誤差を抑えやすいためです。
SciPyのbutter
と
sosfilt
に形式・実装の説明があります。
ziを渡し、zfを次のziにする
最小限のループは次の形です。blocks は、同じ連続信号を時系列順に分割した非空の配列を返すイテレータとします。
import numpy as np
from scipy import signal
sos = signal.butter(6, 40, fs=1000, output="sos")
state = np.zeros((len(sos), 2), dtype=np.float64)
for block in blocks:
block = np.asarray(block, dtype=np.float64)
filtered, state = signal.sosfilt(sos, block, zi=state)
consume(filtered)
blocks と consume は取得系・出力系に合わせて置き換えます。zi を省略すると初期状態はゼロで、戻り値は出力だけです。zi を指定すると (出力, 最終状態) を受け取れます。返された状態を更新し続けることが重要です。
ブロックの長さは一定でなくても構いません。最後が128サンプル未満でも、順序と状態を保てば処理できます。空ブロックは呼び出し前に読み飛ばしてください。
一括処理との一致を数値で確かめる
次は単独で実行できる検証です。1 kHzで2秒間の信号を作り、6次・40 Hzのローパスを適用します。入力はDC成分と8 Hz、180 Hzの正弦波なので乱数による差はありません。
import numpy as np
from scipy import signal
fs = 1000.0
t = np.arange(2000) / fs
x = 1.5 + np.sin(2 * np.pi * 8 * t) + 0.3 * np.sin(2 * np.pi * 180 * t)
sos = signal.butter(6, 40, fs=fs, output="sos")
zi0 = signal.sosfilt_zi(sos) * x[0]
y_whole, z_whole = signal.sosfilt(sos, x, zi=zi0.copy())
block_size = 128
state = zi0.copy()
parts, reset_parts = [], []
for start in range(0, len(x), block_size):
block = x[start:start + block_size]
y_block, state = signal.sosfilt(sos, block, zi=state)
parts.append(y_block)
# 比較用の誤った実装: 各ブロックでゼロ状態に戻す
reset_parts.append(signal.sosfilt(sos, block))
y_stream = np.concatenate(parts)
y_reset = np.concatenate(reset_parts)
np.testing.assert_allclose(y_stream, y_whole, rtol=1e-13, atol=1e-13)
np.testing.assert_allclose(state, z_whole, rtol=1e-13, atol=1e-13)
print(f"max streaming error: {np.max(np.abs(y_stream - y_whole)):.3e}")
error = np.max(np.abs(y_reset[block_size:] - y_whole[block_size:]))
print(f"max reset error after first block: {error:.6f}")
NumPy 2.5.4・SciPy 1.18.1で確認した実行結果は下記のとおりです。最初のブロックを除いても、状態を捨てる実装には大きな誤差が残ります。
max streaming error: 0.000e+00
max reset error after first block: 2.499462
同じ係数・データ型・初期状態を使う今回の検証では出力と最終状態が一致しました。別の実装・演算順序では丸め差があり得るため、検証は許容誤差付きで行っています。比較対象を既定のゼロ状態で一括処理すると、初期条件が違うので一致しません。
破線はブロック境界です。状態を引き継ぐ出力は一括処理と重なります。ゼロに戻す出力は各境界で立ち上がり直します。波形によっては段差、音声ならクリック音として現れます。図を再生成する完全なコードは 検証スクリプト にあります。
sosfilt_ziで初期化するときの仮定
sosfilt_zi(sos) * x[0] は、開始前から入力が最初の値で一定だったとみなす初期化です。
sosfilt_ziの公式説明
では、単位ステップに対する定常状態を求めています。
今回のローパスはDCゲインが1なので、初期値1.5の一定入力に対する起動時の立ち上がりを避けられます。しかし、未知の過去の正弦波の状態を復元する方法ではありません。正弦波や急変が開始直前に存在していれば、実際の履歴との差による過渡応答は残ります。
| 初期化 | 想定している履歴 | 適した場面 |
|---|---|---|
| ゼロ状態 | 過去の入力と状態がゼロ | 信号がゼロから実際に開始する |
sosfilt_zi(sos) * x[0] | 過去の入力が最初の値で一定 | DCを含む連続計測の開始 |
前回の zf | 過去の処理履歴を継続 | 同じストリームの次のブロック |
| 取得済みの先行データでウォームアップ | その先行データの履歴 | 記録途中からの評価開始 |
ハイパス・バンドパスではDCが通らないため、定常出力は一般に入力値そのものではありません。低域通過の例をそのまま「出力は必ず最初の値になる」と一般化しないでください。
多チャンネルではaxisと状態の形をそろえる
入力を (サンプル数, チャンネル数) にした場合、時間方向は axis=0 です。状態は時間軸の長さを2に置き換え、その先頭にセクション軸を加えるので (K, 2, チャンネル数) になります。
# 前の検証コードに続けて実行
X = np.column_stack([x, 2 * x]) # (2000, 2)
zi = signal.sosfilt_zi(sos)[:, :, None] * X[0][None, None, :]
Y, zf = signal.sosfilt(sos, X, axis=0, zi=zi)
print(zf.shape) # (3, 2, 2)
np.testing.assert_allclose(Y[:, 0], y_whole, rtol=1e-13, atol=1e-13)
np.testing.assert_allclose(Y[:, 1], 2 * y_whole, rtol=1e-13, atol=1e-13)
入力が (チャンネル数, サンプル数) なら axis=-1、状態は (K, チャンネル数, 2) です。チャンネル間で同じ状態を共有すると、別の信号の履歴が混ざります。2チャンネルでは形だけで間違いを検出しにくいため、各列の独立処理と一致することも確認してください。
係数・入力・状態をまず float64 でそろえ、想定する最大振幅でも有限値が返るか検証します。float32 や固定小数点を使う場合は、高次数・狭帯域・長時間のデータで精度を評価してください。SOSにしただけで任意の量子化や飽和に耐えられるわけではありません。
sosfiltfiltに置き換えると同じ特性にはならない
オフラインの sosfiltfilt は前後方向に処理します。実係数フィルタの理想的な内部区間では、合成周波数応答は
\[ H_{\mathrm{fb}}(e^{j\omega}) = H(e^{j\omega})H(e^{-j\omega}) = |H(e^{j\omega})|^2 \]となり、位相が相殺されます。振幅ゲイン自体が一方向処理の二乗になるので、遮断周波数で一方向が約−3.01 dBなら前後方向は約−6.02 dBです。これは電力ゲインと振幅ゲインを取り違えた話ではありません。
# 前の検証コードに続けて実行
_, h = signal.sosfreqz(sos, worN=[40], fs=fs)
print(f"one pass: {20 * np.log10(abs(h[0])):.4f} dB")
print(f"forward-backward: {20 * np.log10(abs(h[0]) ** 2):.4f} dB")
# one pass: -3.0103 dB
# forward-backward: -6.0206 dB
この式は端点処理を無視できる周波数応答の説明です。有限長の記録ではパディング・初期化の影響が端点に残ります。128サンプルずつ sosfiltfilt を適用しても、一括のゼロ位相処理にはなりません。また、同じ6次フィルタを2回かけた特性は、12次バターワースを一方向にかけた特性とも一般には異なります。
filtfiltの説明
も参照してください。
ブロック待ち時間と群遅延を分けて考える
128サンプルを1 kHzで集めてから処理するなら、ブロックの時間幅は128 msです。先頭のサンプルは、最後のサンプルが届くまで約127 ms待ちます。これに取得バッファ、スケジューリング、計算、出力の待ち時間が加わります。
一方、IIRの群遅延は位相特性の周波数微分で、周波数によって異なります。
\[ \tau_g(\omega) = -\frac{d\arg H(e^{j\omega})}{d\omega} \]ここで \(\omega\) はrad/sample、\(\tau_g\) はサンプル単位です。秒へはサンプリング周波数で割ります。ブロックを小さくすれば取得待ちは短くできますが、フィルタの係数が同じなら群遅延は変わりません。位相遅れや過渡応答を、固定の128 msだけで説明しないことが重要です。
SciPyのコード例は処理結果の検証用であり、OS上のコールバック期限を保証するものではありません。音声用途ではコールバック内の配列確保やファイル書き込みを避け、実機で処理時間の最大値も測定してください。
実装で確認する項目
- 同じ連続ストリームでは、
zfを次のziに渡す。 - チャンネルごとの状態を独立に保ち、時間軸を明示する。
- 一括処理との比較では、同じ初期条件を使う。
- 欠落サンプルは「なかったこと」にせず、時刻・欠落数と補完方針を管理する。
- 再接続や係数変更の際は初期化方針を決める。別の係数の状態をそのまま流用すると、過渡変化が起きることがある。
よくある質問
sosfiltでブロック処理するときに必要なものは？
前のブロックで返された最終状態zfを、次のブロックの初期状態ziとして渡します。係数と初期状態が同じなら、連続した信号の分割処理は一括処理と一致します。
sosfiltfiltはリアルタイム処理に使えますか？
通常の前後方向処理は未来のサンプルを必要とするため、遅延なしの因果的ストリーミングには使えません。ブロックごとに適用すると端点処理も変わり、全体の一括処理とは一致しません。