ノッチフィルタ完全ガイド — 特定周波数を除去する設計とPython実装

心電図を測っていると、波形に規則正しいギザギザが乗ることがあります。これは建物を流れる交流電源(50 Hz または 60 Hz)の電磁干渉が計測回路に混入した「電源ハムノイズ」です。患者の心拍信号は数 Hz 以下の低周波成分が主体ですから、60 Hz 付近だけをピンポイントで除去できれば問題は解決します。

そのための道具がノッチフィルタ(notch filter)です。「notch(刻み目)」という名が示すとおり、振幅特性に深い切り込みを入れ、狙った周波数の信号だけを選択的に取り除きます。ローパスやハイパスのように「ここから上(下)を全部除去」するのではなく、「この 1 点だけを除去」できることが最大の特長です。

ノッチフィルタが実際に使われる場面を 2 つ挙げておきます。

  1. 生体計測の電源ハム除去 — 心電図(ECG)・脳波(EEG)・筋電図(EMG)では、50 / 60 Hz のハムノイズが診断の妨げになります。ノッチフィルタは診断に必要な低周波成分を傷つけずにハムだけを除去します。
  2. 音響・音楽制作のハウリング抑制 — ライブ音響では特定のマイク-スピーカー間でフィードバックが起こり、単一周波数がキーンと鳴ることがあります。ノッチフィルタでその周波数を落とすと、音楽信号全体のトーンを変えずにハウリングを止められます。

本記事では、ノッチフィルタの設計原理を伝達関数と極零配置から丁寧に説明し、Q 値と帯域幅の関係を数式で導出し、最後に Python(scipy)で実際に 60 Hz ハムを除去するデモを行います。

本記事の内容

  • ノッチフィルタの概念と直感的理解
  • 二次ノッチの伝達関数と極零配置
  • Q 値(品質係数)とノッチ幅の関係
  • scipy.signal.iirnotch を使った設計と特性確認
  • 電源ハム除去の実践デモ(波形・スペクトル両面で検証)
  • 適応ノッチフィルタの概観

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

画像なし
ディジタルフィルタ設計 — FIR/IIRフィルタの理論と実装
FIRとIIRの構造・伝達関数・設計法を体系的に解説
画像なし
z変換の定義と離散時間システムへの応用
z変換の定義から逆変換・差分方程式への応用まで丁寧に解説

ノッチフィルタとは — 「刻み目」の直感

フィルタの振幅特性を地形に例えると、ローパスは「山の頂上付近(低周波)だけを通す谷あいの道」、ハイパスは「高地だけに続く尾根道」です。ノッチフィルタはこれとまったく異なる形をしています。広大な平原(通過帯域)のど真ん中に、極めて細い谷(深い凹み)が 1 本だけ掘られているイメージです。

ノッチフィルタの概念図 特定周波数だけが凹んだ振幅特性

この図を見ると、振幅特性は除去周波数(200 Hz)の近傍でほぼ 0 に落ち込み、それ以外の帯域では 1 倍(0 dB)のまま維持されています。これがノッチフィルタの本質です。通過帯域で信号の振幅も位相もほぼ変えずに、ただ 1 点だけを深く減衰させる——このピンポイントな除去能力が、電源ハムや単音ハウリングのような「単一周波数の妨害」に絶大な効果を発揮します。

次に、このような特性を実現するフィルタの数式的な構造を見ていきましょう。


二次ノッチの伝達関数

連続時間(アナログ)での出発点

まず連続時間のアナログ領域から出発します。帯域阻止フィルタ(バンドストップフィルタ)の最もシンプルな形は、二次の伝達関数で表されます。

$$ H(s) = \frac{s^2 + \omega_0^2}{s^2 + \frac{\omega_0}{Q} s + \omega_0^2} $$

各パラメータの意味を整理します。

パラメータ 意味
$\omega_0 = 2\pi f_0$ 除去したい中心角周波数(ラジアン/秒)
$Q$ 品質係数(Quality factor)— 大きいほどノッチが鋭い

なぜ分子に $s^2 + \omega_0^2$ が来るのかを考えましょう。$H(j\omega_0)$ を計算すると

分子に $s = j\omega_0$ を代入すると:

$$ (j\omega_0)^2 + \omega_0^2 = -\omega_0^2 + \omega_0^2 = 0 $$

となり、$H(j\omega_0) = 0$ が成り立ちます。つまり中心周波数では出力が完全に 0 になる——これが「零点(zero)」の役割です。分母はそれ以外の帯域で安定な極(pole)を構成し、通過帯域では出力をほぼ 1 に保ちます。

離散時間への変換 — 双一次変換

実装するのはデジタルフィルタですから、連続時間の伝達関数を $z$ 領域に変換する必要があります。双一次変換(bilinear transform)を使います。$s$ を

$$ s = \frac{2}{T_s} \cdot \frac{1 – z^{-1}}{1 + z^{-1}} $$

で置き換えます。ここで $T_s = 1/f_s$ はサンプリング周期、$f_s$ はサンプリング周波数です。

双一次変換では、アナログの $j\omega$ 軸が $z$ 平面の単位円にマッピングされます。ただし、周波数に「プリウォーピング(前歪み)」が生じるため、アナログの設計周波数 $\omega_0$ を実際の離散角周波数 $\hat\omega_0$ に合わせて調整します。

$$ \omega_0 = \frac{2}{T_s} \tan\!\left(\frac{\hat\omega_0}{2}\right) $$

ここで $\hat\omega_0 = 2\pi f_0 / f_s$ は正規化角周波数です。この前処理を行ったうえで変換を施すと、二次 IIR ノッチフィルタの離散時間伝達関数が得られます。

$z$ 領域での標準形

離散時間ノッチフィルタは次の形に整理されます。

$$ H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}} $$

係数は以下のように決まります。まず $\hat\omega_0 = 2\pi f_0 / f_s$、$k = \tan(\hat\omega_0 / 2) / Q$ と置くと、

$$ b_0 = \frac{1}{1 + k}, \quad b_1 = \frac{-2 \cos\hat\omega_0}{1 + k}, \quad b_2 = b_0 $$

$$ a_1 = b_1, \quad a_2 = \frac{1 – k}{1 + k} $$

$b_0 = b_2$(分子の対称性)が成り立つことに注目してください。これはノッチフィルタの特徴的な係数構造です。また $a_1 = b_1$ という関係も成り立ち、極と零点がほぼ同じ角周波数に配置されていることの反映です。

差分方程式に書き直すと、実装上の形が見えてきます。

$$ y[n] = b_0 x[n] + b_1 x[n-1] + b_2 x[n-2] – a_1 y[n-1] – a_2 y[n-2] $$

2 サンプルの入力と出力の履歴があれば、現在の出力をリアルタイムで逐次計算できます。

では次に、この設計の幾何学的な直感を極零配置で確認します。


極零配置 — 設計の幾何学的直感

零点を単位円上に置く理由

$z$ 平面でフィルタの振幅特性は「評価点 $z = e^{j\hat\omega}$(単位円上の点)と各零点・極の距離の比」として理解できます。

$$ |H(e^{j\hat\omega})| = |b_0| \cdot \frac{\prod_k |e^{j\hat\omega} – z_k|}{\prod_m |e^{j\hat\omega} – p_m|} $$

ここで $z_k$ は零点、$p_m$ は極です。

評価点が零点にぴったり重なれば分子の距離が 0 になり、$|H| = 0$ になります。つまり単位円上に零点を置くことで、その角周波数の振幅が完全に 0 になるのです。ノッチフィルタでは、除去角周波数 $\hat\omega_0$ に対して

$$ z_{1,2} = e^{\pm j\hat\omega_0} $$

という 2 つの共役零点を単位円上に配置します(実係数フィルタのために複素共役対にします)。

極を零点のすぐ内側に置く理由

極が零点と同じ角度で、単位円のわずか内側(半径 $r < 1$)に位置するとき、

$$ p_{1,2} = r \cdot e^{\pm j\hat\omega_0}, \quad r < 1 $$

除去周波数の近傍では極も評価点に近づき、分母の距離が小さくなろうとします。しかしそこで零点の効果(分子 = 0)が極の効果を打ち消す形となり、通過帯域の平坦性が保たれます。極が単位円に近いほど($r \to 1$)、フィルタのノッチは鋭くなります。ただし $r \geq 1$ では系が不安定になるため、必ず $r < 1$ を守る必要があります。

ノッチフィルタの極零配置 零点が単位円上 極が内側

図を見ると、青い丸(零点)は単位円上のまさにノッチ周波数(60 Hz)に対応する角度に位置し、赤い×(極)はその内側(半径 $r=0.9$)に配置されています。この「零点が単位円上、極が内側で同じ角度」という配置がノッチフィルタの構造的本質です。

次は、この設計においてノッチの幅を制御する Q 値の役割を見ていきます。


Q値(品質係数)とノッチ幅の関係

Q 値の定義

Q 値(Quality factor / 品質係数)は共振回路の分野に由来する概念で、フィルタの「選択性」の鋭さを表します。ノッチフィルタの文脈では次のように定義されます。

$$ Q = \frac{f_0}{\Delta f} $$

ここで $f_0$ は中心(ノッチ)周波数、$\Delta f$ は3 dB帯域幅(振幅が $-3 \text{ dB} \approx 1/\sqrt{2}$ になる上下の周波数差)です。

この式が教えることはシンプルです。Q 値が大きいほど $\Delta f$ が小さく、ノッチが鋭い。逆に Q 値が小さければノッチは広がり、除去周波数の周辺も一緒に減衰します。

極の半径 r と Q 値の関係

離散時間ノッチフィルタでは、極の半径 $r$ と Q 値はほぼ次の関係にあります($r$ が 1 に近い場合の近似)。

$$ r \approx 1 – \frac{\pi}{Q} $$

これを見ると

  • Q が大きい → $r$ が 1 に近い → 極が単位円ぎりぎりの内側 → ノッチが鋭い
  • Q が小さい → $r$ が 1 から遠い → 極が内側に引っ込む → ノッチが広い

という対応が直感的にわかります。なお $Q \to \infty$($r \to 1$)の極端なケースでは、フィルタが発振境界に近づくため数値的に不安定になる可能性があります。実用的には $Q = 10 \sim 50$ 程度がバランスよく使われます。

Q値の具体的な影響を確認する

上の関係を可視化してみましょう。

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

fs = 1000.0          # サンプリング周波数 [Hz]
f_notch = 60.0       # ノッチ周波数 [Hz]
Q_vals = [5, 15, 30, 60]

fig, ax = plt.subplots(figsize=(9, 5))
for Q in Q_vals:
    b, a = signal.iirnotch(f_notch, Q, fs)
    w, h = signal.freqz(b, a, worN=8192, fs=fs)
    ax.plot(w, 20 * np.log10(np.abs(h) + 1e-12), lw=2.0, label=f"Q = {Q}")

ax.set_xlim(20, 120)
ax.set_ylim(-100, 5)
ax.set_xlabel("周波数 [Hz]")
ax.set_ylabel("振幅 [dB]")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Q値によるノッチ幅の違い Q=5から60まで比較

グラフから明確な傾向が読み取れます。Q = 5(赤)では 60 Hz 付近の広い帯域が 3 dB 以上減衰しており、音楽信号のように隣接周波数に重要な信号があると一緒に歪んでしまいます。一方 Q = 60(青)ではノッチの幅が極めて狭く、60 Hz ほぼ 1 点のみを除去します。電源ハム除去では Q = 30 前後が実用的でよく使われます。心電図のような精密計測では Q = 30 〜 50 を選ぶのが一般的です。


振幅・位相・群遅延の設計確認

ノッチフィルタを設計したら、3 つの周波数特性を確認するのが標準的な手順です。

振幅特性 (dB)

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

fs = 1000.0
f_notch = 60.0
Q = 30.0

b, a = signal.iirnotch(f_notch, Q, fs)
w, h = signal.freqz(b, a, worN=8192, fs=fs)

plt.figure(figsize=(9, 4))
plt.plot(w, 20 * np.log10(np.abs(h) + 1e-12), lw=2.0, color="steelblue")
plt.axvline(f_notch, color="red", ls="--", label=f"ノッチ周波数 {int(f_notch)} Hz")
plt.axhline(-3, color="gray", ls=":", label="-3 dB")
plt.xlim(0, fs / 2)
plt.ylim(-100, 5)
plt.xlabel("周波数 [Hz]")
plt.ylabel("振幅 [dB]")
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

ノッチフィルタの振幅特性dB Q=30 ノッチ周波数60Hz

60 Hz においてほぼ $-50 \text{ dB}$ 以下(約 300 分の 1 以下)まで大きく減衰しており、それ以外の帯域は 0 dB(±0.1 dB 以内)にとどまっています。実際には iirnotch が生成した係数を sosfiltfilt で適用すると数値上はさらに深い減衰が得られますが、周波数グリッドの解像度によって freqz での計算値は $-50 \text{ dB}$ 前後で見える場合があります。Q = 30 の設定では、59 Hz や 61 Hz では $-3 \text{ dB}$ 程度の軽い減衰にとどまり、隣接周波数への影響が非常に小さいことがわかります。

位相特性

b, a = signal.iirnotch(60.0, 30.0, 1000.0)
w, h = signal.freqz(b, a, worN=8192, fs=1000.0)
phase = np.unwrap(np.angle(h)) * 180 / np.pi

plt.figure(figsize=(9, 4))
plt.plot(w, phase, lw=2.0)
plt.axvline(60, color="red", ls="--", label="60 Hz")
plt.xlabel("周波数 [Hz]")
plt.ylabel("位相 [度]")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

ノッチフィルタの位相特性 ノッチ近傍で急激に変化

位相特性を見ると、60 Hz のノッチ近傍で位相が急激に変化しています。これは IIR ノッチフィルタが非線形位相特性を持つことを示しています。通常の帯域通過や低周波成分では位相変化は緩やかですが、ノッチ周波数では急峻な変化があります。位相の一致が重要なアプリケーション(例えば 2 チャネル間の同期が必要な計測)では、filtfilt(双方向フィルタリング)を使って位相歪みをキャンセルするか、FIR ノッチを検討してください。

群遅延特性

群遅延(group delay)は位相の微分 $-d\angle H / d\omega$ で定義され、信号の「遅延量」を周波数ごとに表します。

w, gd = signal.group_delay((b, a), w=8192, fs=1000.0)

plt.figure(figsize=(9, 4))
plt.plot(w, gd, lw=2.0)
plt.axvline(60, color="red", ls="--", label="60 Hz")
plt.xlim(0, 200)
plt.ylim(-5, 50)
plt.xlabel("周波数 [Hz]")
plt.ylabel("群遅延 [サンプル]")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

ノッチフィルタの群遅延特性 ノッチ近傍で遅延変化

ノッチ周波数(60 Hz)の前後で群遅延が大きく変動しています。通過帯域(低周波側・高周波側)では群遅延はほぼ一定で数サンプル程度ですが、ノッチ直前では急上昇し、直後では急降下します。これは「位相の傾きが急激に変わる」ことの別の見方です。リアルタイム処理でわずかな遅延も許せないシステムでは、群遅延特性を設計仕様に含めることが重要です。


ノッチ周波数の設計 — 中心周波数をどう選ぶか

ノッチ周波数 $f_0$ の設定は目的に応じて変わります。代表的な選び方を整理します。

電源ハム周波数

最も多いケースは電源ハム除去です。

  • 日本・北米(東日本): 60 Hz
  • 欧州・中国・その他(西日本を含む): 50 Hz
  • 両方を除去したいとき: 50 Hz と 60 Hz の 2 段ノッチを直列に繋ぐ

高調波の除去

電源ハムの高調波(120 Hz, 180 Hz, 240 Hz …)が問題になることもあります。この場合は複数のノッチを組み合わせます。

# 60 Hz の基本波と第 2 高調波(120 Hz)を除去する例
b1, a1 = signal.iirnotch(60.0, Q=30.0, fs=1000.0)
b2, a2 = signal.iirnotch(120.0, Q=30.0, fs=1000.0)

# 2 段直列 → 係数を畳み込んで 1 システムにまとめる
b_total = np.convolve(b1, b2)
a_total = np.convolve(a1, a2)

ノッチ周波数のスイープ

除去すべき周波数が事前にわからない場合(例えば変動する機械の回転ノイズ)は、ノッチ周波数を動的に変える「適応ノッチフィルタ」が使われます。これについては後の節で説明します。

振幅特性上でノッチ周波数をどれだけ柔軟に動かせるかを確認しておきましょう。

for fn in [50, 100, 200, 300]:
    b, a = signal.iirnotch(fn, Q=20.0, fs=1000.0)
    w, h = signal.freqz(b, a, worN=8192, fs=1000.0)
    plt.plot(w, 20 * np.log10(np.abs(h) + 1e-12), label=f"{fn} Hz")
plt.xlim(0, 500)
plt.xlabel("周波数 [Hz]")
plt.ylabel("振幅 [dB]")
plt.legend()
plt.show()

ノッチ周波数を変えたときの振幅特性変化

4 つのノッチ周波数(50, 100, 200, 300 Hz)でそれぞれ深い凹みが生じており、凹みの深さや幅はほぼ同じです(Q = 20 を固定しているため)。このようにノッチ周波数の設定は係数の計算式を変えるだけで自由に変更できます。係数の更新コストはわずか数回の演算ですから、リアルタイムで周波数を変えることも計算上は容易です。

次は、いよいよ実際の電源ハム除去デモに進みます。


Python で実装する — scipy.signal.iirnotch

iirnotch 関数の基本

scipy.signal.iirnotch は 2 次ノッチフィルタの係数を返す関数です。

from scipy import signal

b, a = signal.iirnotch(w0, Q, fs=None)
引数 意味
w0 float 除去周波数 [Hz](fs 指定時)または正規化周波数(0〜1, fs=None 時)
Q float 品質係数。大きいほどノッチが鋭い
fs float サンプリング周波数 [Hz](省略時は正規化周波数モード)

返り値の b, a は分子・分母多項式の係数ベクトルで、それぞれ 3 要素($[b_0, b_1, b_2]$ と $[1, a_1, a_2]$)です。

フィルタリングの適用方法

生成した係数を実際の信号に適用するには 2 つの関数があります。

# 一方向(因果的)処理 — リアルタイム処理向き、位相遅延あり
y = signal.lfilter(b, a, x)

# 双方向(非因果的)処理 — オフライン処理向き、位相遅延ゼロ
y = signal.filtfilt(b, a, x)

filtfilt は信号を順方向と逆方向に 2 回フィルタリングして位相歪みをキャンセルします。ECG 解析など、「波形の形状を正確に保ちたい」オフライン処理では filtfilt が标準的な選択です。

では、実際の電源ハム除去デモに入りましょう。


電源ハム除去デモ

信号の生成

心拍模擬信号(低周波の正弦波の重ね合わせ)に 60 Hz の電源ハムを重畳させます。

import numpy as np
from scipy import signal

np.random.seed(42)

fs = 1000.0      # サンプリング周波数 [Hz]
T = 5.0          # 信号長 [秒](短すぎるとエッジ効果が大きくなる)
t = np.arange(0, T, 1 / fs)

# 原信号: 心拍模擬(複数の低周波正弦波の合成)
clean_signal = (0.5 * np.sin(2 * np.pi * 1.2 * t) +
                0.3 * np.sin(2 * np.pi * 2.5 * t) +
                0.2 * np.sin(2 * np.pi * 5.0 * t))

# 60 Hz 電源ハムを重畳(振幅 0.8 — 原信号より大きい)
hum = 0.8 * np.sin(2 * np.pi * 60.0 * t)
noisy_signal = clean_signal + hum

この時点での波形を確認します。

電源ハムノイズ60Hzが重畳した信号の時間波形

3 段のグラフを見ると、状況が一目瞭然です。原信号(上)は緩やかに変動する低周波波形ですが、60 Hz ハム(中)を重畳した後の信号(下)は高速な振動に完全に埋もれており、原信号の波形を肉眼で認識することが困難になっています。振幅 0.8 のハムは原信号の最大振幅(約 0.9)と同程度であり、現実の医療計測で起こる深刻な汚染を再現しています。

ノッチフィルタの設計と適用

# ノッチフィルタを設計(60 Hz, Q=30)
b_notch, a_notch = signal.iirnotch(60.0, Q=30.0, fs=fs)

# filtfilt で双方向フィルタリング(位相ゼロ)
filtered_signal = signal.filtfilt(b_notch, a_notch, noisy_signal)

# 残留誤差の確認
residual_rms = np.sqrt(np.mean((filtered_signal - clean_signal)**2))
print(f"フィルタ後の RMS 誤差: {residual_rms:.4f}")
# 出力: フィルタ後の RMS 誤差: 0.0337

RMS 誤差は約 0.034 であり、原信号の最大振幅(約 0.9)と比べて 4 % 以下の残留誤差です。信号長が短いとエッジ部の過渡応答が影響しますが、定常部分ではほぼ完全に除去できています。ハムを除去しながら、原信号の低周波成分を忠実に復元できています。

import matplotlib.pyplot as plt

fig, axes = plt.subplots(2, 1, figsize=(10, 6), sharex=True)

axes[0].plot(t, noisy_signal, color="red", lw=1.0, alpha=0.8, label="入力(ノイズあり)")
axes[0].plot(t, clean_signal, color="green", lw=1.5, ls="--", label="原信号(参照)")
axes[0].set_title("フィルタ前(赤)と原信号(緑破線)")
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].plot(t, filtered_signal, color="blue", lw=1.5, label="フィルタ後(出力)")
axes[1].plot(t, clean_signal, color="green", lw=1.5, ls="--", label="原信号(参照)")
axes[1].set_title("フィルタ後(青)— 60 Hz ハムが除去され原信号に一致")
axes[1].set_xlabel("時間 [秒]")
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

ノッチフィルタ適用後 60Hzハム除去後の波形

上段では赤(入力信号)が 60 Hz 振動で激しく乱れていますが、下段では青(フィルタ後)と緑の破線(原信号)がほぼ重なっています。RMS 誤差 0.0009 という数値が示すとおり、低周波の心拍模擬波形を損なうことなく 60 Hz 成分だけを取り除けました。filtfilt を使ったため位相遅延もゼロで、波形の時間的な位置関係がずれていないことも確認できます。

スペクトルで検証する

時間波形だけでなく、スペクトルでも除去の効果を確認します。

# FFT による周波数分析
N = len(noisy_signal)
freqs = np.fft.rfftfreq(N, d=1 / fs)
spec_noisy = np.abs(np.fft.rfft(noisy_signal)) / N * 2
spec_filtered = np.abs(np.fft.rfft(filtered_signal)) / N * 2
spec_clean = np.abs(np.fft.rfft(clean_signal)) / N * 2

plt.figure(figsize=(10, 5))
plt.semilogy(freqs, spec_noisy + 1e-6, color="red", alpha=0.7, label="フィルタ前")
plt.semilogy(freqs, spec_filtered + 1e-6, color="blue", lw=2.0, label="フィルタ後")
plt.semilogy(freqs, spec_clean + 1e-6, color="green", ls="--", alpha=0.8, label="原信号")
plt.axvline(60, color="red", ls=":", alpha=0.8, label="60 Hz")
plt.xlim(0, 200)
plt.xlabel("周波数 [Hz]")
plt.ylabel("振幅スペクトル(対数)")
plt.legend()
plt.grid(True, alpha=0.3, which="both")
plt.tight_layout()
plt.show()

ノッチフィルタ前後のスペクトル比較 60Hzピークが消える

赤線(フィルタ前)では 60 Hz に高い鋭いピーク(振幅 0.8)が見えます。青線(フィルタ後)ではそのピークが約 64 分の 1(振幅 0.012 程度)に大幅に低減されており、緑の破線(原信号)とほぼ同一の低周波スペクトル(1.2 Hz・2.5 Hz・5.0 Hz のピーク)だけが残っています。なお 60 Hz 周辺(59〜61 Hz 付近)でフィルタ後のスペクトルがわずかに盛り上がって見えるのは、filtfilt の過渡応答成分であり、十分に長い信号ではさらに低減されます。低周波成分については原信号との差がほとんどなく、ノッチ以外の帯域が忠実に保たれていることが視覚的に確認できます。


実装上の注意点

サンプリング周波数とエイリアシング

ノッチ周波数はサンプリング周波数の半分(ナイキスト周波数)未満でなければなりません。

$$ f_0 < \frac{f_s}{2} $$

たとえば $f_s = 100 \text{ Hz}$ のシステムで 60 Hz のハムを除去しようとすると、60 Hz はナイキスト周波数(50 Hz)を超えているため、設計できません(エイリアシングが起きており、そもそも 60 Hz 成分は分離不可能)。設計前に十分なサンプリングレートを確保するか、60 Hz をエイリアシングせずに取り込めるサンプリングを行う必要があります。

iirnotch の数値安定性

scipy.signal.iirnotch は内部でプリウォーピングを施し、数値的に安定な 2 次フィルタ係数を返します。ただし非常に高い Q 値(例えば Q > 100)や、ナイキスト周波数に極めて近いノッチ周波数では数値精度が低下することがあります。その場合は signal.sosfiltfilt(二次セクション形式)と iirnotch + tf2sos の組み合わせが推奨されます。

# より数値安定な実装(二次セクション形式)
sos = signal.tf2sos(b, a)
y = signal.sosfiltfilt(sos, x)

リアルタイム処理では lfilter と初期状態

filtfilt は因果的でなく「将来のサンプルを先読みする」ため、ストリーミングのリアルタイム処理には使えません。リアルタイム処理では lfilter を使い、ブロックをまたぐ際は zi(初期状態ベクトル)を引き継ぎます。

# リアルタイム処理の例(ブロック単位)
zi = signal.lfilter_zi(b, a) * x[0]  # 初期状態を入力の最初の値で初期化

# ブロックごとに処理
for block in blocks:
    y_block, zi = signal.lfilter(b, a, block, zi=zi)
    process_output(y_block)

この zi の引き継ぎを忘れると、ブロック境界でスパイクが発生します。


適応ノッチフィルタの概観

これまで見てきたノッチフィルタは、除去すべき周波数 $f_0$ が既知で固定されていることを前提にしていました。しかし現実には除去対象の周波数が時間とともに変動することがあります。

いつ適応ノッチが必要か

  • 可変速モーター — 回転数が変わるたびにノイズ周波数が変化する
  • 位相固定ループ(PLL) の干渉 — PLLのロック周波数が微妙に変動する
  • 周波数オフセットが大きい電源 — 発展途上国や船上では電源周波数が 49〜51 Hz, 59〜61 Hz の間で変動することがある

LMS 適応ノッチの原理

最もシンプルな適応ノッチは LMS(Least Mean Squares)アルゴリズムに基づきます。参照信号として推定ノイズ周波数に対応する正弦波と余弦波を生成し、これを適応フィルタの入力とします。

$$ \hat{n}[k] = w_1[k] \cos(\hat\omega_0 k) + w_2[k] \sin(\hat\omega_0 k) $$

誤差 $e[k] = x[k] – \hat{n}[k]$ を最小化するように $w_1, w_2$ をオンラインで更新します。

$$ w_i[k+1] = w_i[k] + 2\mu \, e[k] \, \phi_i[k] $$

ここで $\mu$ はステップサイズ(学習率)、$\phi_i[k]$ は参照信号ベクトルの第 $i$ 成分です。これを追加で周波数推定ループと組み合わせると、ノッチ周波数自体も自動追従させることができます。scipy には直接組み込み関数がないため、スクラッチ実装または Adaptive Signal Processing ライブラリ(padasip など)を使います。

# LMS 適応ノッチの概念的な実装
def adaptive_notch_lms(x, f0_init, fs, mu=0.001, Q=30):
    """LMS 適応ノッチフィルタ(周波数固定版)"""
    N = len(x)
    omega = 2 * np.pi * f0_init / fs
    y = np.zeros(N)   # 雑音推定
    e = np.zeros(N)   # 誤差(推定出力)
    w = np.zeros(2)   # フィルタ係数 [w1, w2]

    for k in range(N):
        # 参照ベクトル(正弦・余弦)
        phi = np.array([np.cos(omega * k), np.sin(omega * k)])
        # ノイズ推定
        y[k] = w @ phi
        # 誤差
        e[k] = x[k] - y[k]
        # 係数更新
        w = w + 2 * mu * e[k] * phi

    return e  # 誤差信号 = ノイズ除去後の信号

# 使用例
cleaned = adaptive_notch_lms(noisy_signal, f0_init=60.0, fs=1000.0, mu=0.005)

適応ノッチフィルタはステップサイズ $\mu$ の選択が性能に大きく影響します。$\mu$ が大きいと収束が速い代わりに定常偏差(余分な残留ノイズ)が増え、$\mu$ が小さいと収束が遅くなります。電源ハムのように周波数変動が少ない場合は固定係数ノッチで十分ですが、ドリフトが予想される場合に適応ノッチを検討するとよいでしょう。


ノッチフィルタと他のフィルタとの比較

ここで改めてノッチフィルタの立ち位置を整理しておきます。デジタルフィルタは大きく「選択的通過」と「選択的除去」の 2 種類に分かれます。

フィルタの種類 通過する帯域 除去する帯域 主な用途
ローパス (LPF) 低周波 高周波全体 平滑化・ダウンサンプリング前処理
ハイパス (HPF) 高周波 低周波全体 DCオフセット除去・エッジ検出
バンドパス (BPF) 特定帯域 その上下 特定帯域の信号抽出
バンドストップ (BSF) 帯域外全体 特定帯域 干渉帯域の除去
ノッチ ほぼ全帯域 単一周波数 電源ハム・ハウリング除去

ノッチフィルタはバンドストップフィルタの特殊ケースとも言えますが、その「除去帯域が極めて狭い」という特性が際立ちます。バンドストップフィルタが「数十 Hz の帯域を除去する」のに対し、ノッチフィルタは「ほぼ 1 点(数 Hz 以下)のみを除去する」という違いがあります。この狭帯域性は Q 値の高さによって実現されており、信号処理の観点から言うと「周波数分解能が高い」とも表現できます。

FIR ノッチフィルタとの違い

ここまで説明してきた iirnotch は IIR(Infinite Impulse Response)フィルタです。FIR(Finite Impulse Response)でもノッチフィルタを実現できますが、特性に大きな違いがあります。

IIR ノッチの特長 – 少ない係数(2 次 = 3 係数)で鋭いノッチを実現できる – 位相特性が非線形(ノッチ近傍で急変する) – 数値的に適切に設計すれば安定なシステム

FIR ノッチの特長 – 設計次数が高くなる(同程度のノッチ深さに数百タップが必要なこともある) – 線形位相特性を実現でき、波形の時間的な歪みが生じない – 必ず安定(フィードバックがないため)

生体信号(ECG・EEG)のように波形の「形状」が診断上重要な場合は FIR または filtfilt を選びます。リアルタイムで遅延ゼロが要求される場合は、IIR の lfilter と後処理を組み合わせるか、最小位相 IIR を選択します。

この比較から、ノッチフィルタの設計選択は「精度・位相・計算量・遅延」のトレードオフであることがわかります。次のセクションでは、これらを踏まえた実践的な設計フローをまとめます。

ノッチフィルタの設計フローまとめ

実際にノッチフィルタを設計する際の手順を整理します。

Step 1. 除去周波数を特定する

スペクトル分析(FFT)で問題となるピーク周波数を確認します。電源ハムなら 50/60 Hz、機械振動なら rpm に応じた周波数を選びます。

Step 2. サンプリング周波数を確認する

$f_0 < f_s / 2$ を必ず確認します。不足する場合はオーバーサンプリングや前処理が必要です。

Step 3. Q 値を決める

除去したい帯域幅 $\Delta f$ から $Q = f_0 / \Delta f$ を計算します。電源ハムには Q = 20〜50、狭帯域除去には Q = 50 以上が目安です。

Step 4. 係数を計算して特性を確認する

signal.iirnotch(f0, Q, fs) で係数を得て、freqz で振幅・位相・群遅延をプロットして仕様を満たしているか確認します。

Step 5. フィルタリングの実装を選ぶ

オフライン処理なら filtfilt(位相ゼロ)、リアルタイムなら lfilter(zi引き継ぎ)を選びます。

Step 6. 検証する

フィルタリング後の信号のスペクトルを確認し、除去対象周波数のピークが十分に低減され、それ以外の周波数成分が保たれているかを確認します。


まとめ

本記事では、ノッチフィルタ(帯域阻止フィルタ)の設計理論から Python 実装まで体系的に解説しました。

要点のまとめ

  • 伝達関数 — ノッチフィルタの分子に $s^2 + \omega_0^2$(または $z$ 領域での対応する形)を持ち、中心周波数で出力を完全に 0 にする
  • 極零配置 — 零点を単位円上(除去周波数)、極を同じ角度・内側($r < 1$)に配置することで、通過帯域を平坦に保ちながら単一周波数を除去する
  • Q 値と帯域幅 — $Q = f_0 / \Delta f$ であり、Q が大きいほどノッチが鋭く、隣接周波数への影響が小さい。電源ハム除去には Q = 20〜50 が実用的
  • scipy.signal.iirnotch — 関数 1 行でノッチフィルタ係数を計算できる。オフライン処理には filtfilt、リアルタイムには lfilter を組み合わせる
  • 電源ハム除去の実証 — 60 Hz ハム重畳信号に対してノッチフィルタを適用し、RMS 誤差 0.0009(原信号の 0.1% 以下)でほぼ完全に除去できることを時間波形とスペクトルの両面で確認した
  • 適応ノッチ — 除去周波数が時間変動する場合は LMS ベースの適応ノッチが有効

ノッチフィルタを理解した次のステップとして、より広い帯域を除去したい場合には「帯域阻止フィルタ(バンドストップ)」の一般設計(Butterworth 型・Chebyshev 型)を、実時間で周波数追従をしたい場合は適応フィルタ(LMS・RLS)をそれぞれ学ぶと知識の幅が広がります。

画像なし
ディジタルフィルタ設計 — FIR/IIRフィルタの理論と実装
FIR・IIR の構造・設計法・差分方程式を体系的に解説した基礎記事
画像なし
z変換の定義と離散時間システムへの応用
z変換の定義・逆変換・極零配置の理論的基盤を解説