🌠
符号付きの周波数を通す IIR Filter

2025/08/15に公開
 話すことIIR フィルタは周波数の符号に関係なく処理がされます。本記事では符号を考慮し、すなわち周波数の正か負を指定しフィルタをかけます。(\cos(\omega t) = \frac{1}{2}\left(e^{i\omega t} + e^{-i\omega t} \right) に正の周波数を通すフィルタをかけると \frac{1}{2}e^{i\omega t} を得ます)
イメージはヒルベルト変換をして得る解析信号に近いです。解析信号は正の周波数のみを持ちますが、今回は複素信号を対象にしてフィルタをかけます。

実装として NumPy, SciPy を使います。

 話さないこと今回使用するフィルタは実係数を扱います。複素係数のフィルタを使い符号を意識したフィルタはできるかもしれませんが、今回は扱いません。

 具体例次式で定義される正の周波数と負の周波数が混在した複素信号を例にします。周波数 1.5, 6.0, -3.0 [Hz] を持ちます。
 x(t) = 0.8 \ e^{1.5 \times 2\pi it} + 0.5 \ e^{6.0 \times 2\pi it} + 0.5 \ e^{-3.0 \times 2\pi it}画像は元の信号の周波数ごとの強度と正の周波数を通すフィルタをかけた信号の強度です。青色が元の信号, オレンジ色がフィルタ後の信号の結果です。フィルタをかけることで負の周波数が減衰していることが分かります。

 準備アルゴリズム説明と実装を話す前にヒルベルト変換と正 (または負) の周波数を得るための演算方法を説明します。

 1. ヒルベルト変換記事 を参考にし、ヒルベルト変換後の信号 \hat{x}(t) をフーリエ変換すると以下式になります。\text{sgn(.)} は入力が正の時 1, 負の時に -1 を出力する関数です。
 \hat{X}(f) = -i. \text{sgn}(f) X(f)元の信号を x(t) = e^{i\omega t} とすると、ヒルベルト変換記号 \mathcal{H} を使い以下を得ます。

\begin{align*}
    \mathcal{H}(e^{i\omega t}) &= -i. \text{sgn}(\omega) e^{i\omega t} \\
      &= \begin{cases}
        -i e^{i\omega t} & (\omega \geq 0) \\
        i e^{i\omega t} & (\omega < 0)
      \end{cases} \\
      &= \begin{cases}
        e^{-i\pi/2} e^{i\omega t} & (\omega \geq 0) \\
        e^{i\pi/2} e^{i\omega t} & (\omega < 0)
      \end{cases}
\end{align*}
\omega \geq 0 のとき e^{i\omega t} を -\frac{\pi}{2} 回転し、\omega < 0 のとき \frac{\pi}{2} だけ回転したものになります。

 2. 正 (または負) の周波数を得る演算正の周波数のみを残す下式を定義します。
 y_{+}(t) = \frac{1}{2} x(t) + \frac{i}{2}\mathcal{H}(x(t))x(t) を構成するある角周波数 \omega について考えると、

\begin{align*}
    y_{+}(t) &= \frac{1}{2} e^{i\omega t} + \frac{i}{2} \mathcal{H}(e^{i\omega t}) \\
        &= \begin{cases}
            \frac{1}{2} e^{i\omega t} + \frac{i}{2} (-i e^{i\omega t}) & (\omega \geq 0)\\
            \frac{1}{2} e^{i\omega t} + \frac{i}{2} (i e^{i\omega t}) & (\omega < 0)
        \end{cases} \\
        &= \begin{cases}
            e^{i\omega t} & (\omega \geq 0)\\
            0& (\omega < 0)
        \end{cases} \\
\end{align*}
\omega \geq 0 のときに 0 でない値を持ち、\omega < 0 のときに 0 になります。

前述のヒルベルト変換では \omega \geq 0 のとき e^{i\omega t} を -\frac{\pi}{2} 回転していました。y_{+}(t) の演算では \mathcal{H} の項を \frac{\pi}{2} 回転させるため元の向きに戻ります。そのため \omega \geq 0 のときに 0 でない値を持ちます。\omega < 0 のとき e^{i\omega t} を \frac{\pi}{2} 回転するため y_{+}(t) の演算では \pi 回転したものになり、元の向きと反対になります。このため \omega < 0 のときに 0 になり、負の周波数がなくなります。
負の周波数のみを残す式を下式のように定義します。
 y_{-}(t) = \frac{1}{2} x(t) - \frac{i}{2}\mathcal{H}(x(t))先ほどと同様に考えると以下を得ます。

\begin{align*}
    y_{-}(t) &= \frac{1}{2} e^{i\omega t} - \frac{i}{2} \mathcal{H}(e^{i\omega t}) \\
        &= \begin{cases}
             0 & (\omega \geq 0)\\
             e^{i\omega t} & (\omega < 0)
        \end{cases} \\
\end{align*}
y_{-}(t) の \mathcal{H} の項は -\frac{\pi}{2} だけ回転がかかります。\omega \geq 0 の場合は元の e^{i\omega t} と反対向きになり、\omega < 0 のときは同じ向きになります。

 アルゴリズム・実装
 概要実装ではフィルタを 2つ使い y_{+}, y_{-} を求めます。フィルタを使うことで位相遅れが出るため x(t) をそのまま使う項をフィルタを通したものに置き換えます。(参考論文)
 y_{+}(t) = \frac{1}{2} \mathcal{H}_{0}x(t) + \frac{i}{2}\mathcal{H}_{1}(x(t)) 重要なのはヒルベルト変換で行ったように位相が \frac{\pi}{2} だけ回転することです。(実数係数を持つ) フィルタの性質上、周波数が正と負の時の位相は符号が逆になります。このため \mathcal{H}_{0}, \mathcal{H}_{1} の正の周波数の時に位相の差が \frac{\pi}{2} であれば、負の周波数の時も前述のヒルベルト変換を使ったものを満たします。

 フィルタ作成位相特性のみを変化させるオールパスフィルタを \mathcal{H}_0, \mathcal{H}_1 に使います。オールパスフィルタは次のような形になっています (参考記事)。
H(z)=\prod_{i=1}^{N}\frac{a_i + z^{-1}}{1 + a_i z^{-1}}これらの位相の差を \pi/2 になるように最適化を行います。目的関数は以下になります。k は対象にする範囲内の周波数であり、K はその総数です。
J=\sum_{k=1}^{K} \Big(\mathrm{unwrap}\big(\angle H_0(e^{j\omega_k})-\angle H_1(e^{j\omega_k})\big)-\tfrac{\pi}{2}\Big)^2パラメータ a_i の初期値は 参考記事 より以下を使います。\text{fs} はサンプリング周波数, f_c は対象範囲内の周波数です。
a=(1-\tan(\pi f_c/\text{fs}))/(1+\tan(\pi f_c/\text{fs}))以下のようにフィルタを作成します。_fit_allpass_pair 関数がフィルタ作成関数です。
# 今回使用するライブラリ
import numpy as np
from numpy.typing import NDArray
from scipy import signal, optimize

def _a_from_fc(fc: float, fs: float) -> float:
    """オールパスフィルタの初期値を求める."""
    k = np.tan(np.pi * float(fc) / float(fs))
    a = (1 - k) / (1 + k)
    return float(np.clip(a, -0.999, 0.999))

def _allpass_freqz(a: float, w: NDArray) -> NDArray:
    """オールパスフィルタの各項を求める."""
    z = np.exp(1j * w)
    return (a + z**-1) / (1 + a * z**-1)

def _cascade_freqz(a_list: NDArray, w: NDArray) -> NDArray:
    """オールパスフィルタを求める."""
    H = np.ones_like(w, dtype=complex)
    for a in a_list:
        H *= _allpass_freqz(float(a), w)
    return H

def _fit_allpass_pair(
    fs: float, f1: float, f2: float, nsec: int, wpts: int = 512, init_ratio: float = np.sqrt(2.0)
) -> tuple[NDArray, NDArray]:
     """2つのオールパスフィルタを最適化する.

    Args:
        fs (float): サンプリング周波数
        f1 (float): この範囲で最適化する (下限)
        f2 (float): この範囲で最適化する (上限)
        nsec (int): フィルタの係数の総数
        wpts (int, optional): 周波数範囲内で最適化するサンプル数
        init_ratio (float, optional): 初期値をずらす係数.

    Returns:
        tuple[NDArray, NDArray]: 2つのフィルタ係数
    """
    # 初期値：帯域を等比配置、片方を少しシフト
    f0 = np.geomspace(f1, f2, nsec)
    f1s = np.clip(f0 * init_ratio, f1 * 1.05, f2 / 1.05)
    a0 = np.array([_a_from_fc(fc, fs) for fc in f0], dtype=float)
    a1 = np.array([_a_from_fc(fc, fs) for fc in f1s], dtype=float)

    # 最適化準備
    w = np.linspace(2 * np.pi * f1 / fs, 2 * np.pi * f2 / fs, wpts)
    target = -(np.pi / 2) * np.ones_like(w)
    x0 = np.concatenate([a0, a1])

    # 最適化の目的関数
    def resid(x):
        a0x, a1x = x[:nsec], x[nsec:]
        H0 = _cascade_freqz(a0x, w)
        H1 = _cascade_freqz(a1x, w)
        dphi = np.unwrap(np.angle(H1) - np.angle(H0))
        return dphi - target

    lb = -0.995 * np.ones_like(x0)
    ub = +0.995 * np.ones_like(x0)
    res = optimize.least_squares(resid, x0, bounds=(lb, ub), max_nfev=4000)
    return res.x[:nsec].astype(float), res.x[nsec:].astype(float)

 フィルタをかけるクラス求まったフィルタから正の周波数, 負の周波数を得たい時に以下のような計算ができるクラスを作ります。

\begin{align*}
  y_{+}(t) &= \frac{1}{2} \mathcal{H}_{0}x(t) + \frac{i}{2}\mathcal{H}_{1}(x(t)) \\
  y_{-}(t) &= \frac{1}{2} \mathcal{H}_{0}x(t) - \frac{i}{2}\mathcal{H}_{1}(x(t))
\end{align*}
class OneSidedIIR:
    """IIR 片側通過フィルタ (正 or 負).

    y_pos = 0.5*(H0.x + j.H1.x)
    y_neg = 0.5*(H0.x - j.H1.x)
    """
    def __init__(
            self, fs: float, band: tuple[float, float], nsec: int = 6, mode: str = "pos",
            design: str = "fit", wpts: int = 512) -> None:
        self.fs = fs  # サンプリング周波数
        self.band = band  # この範囲内でフィルタの最適化を行う
        self.nsec = nsec  # フィルタ数
        self.mode = mode  # "pos": 正周波数を抽出, "neg": 負周波数を抽出
        self.design = design  # "fit": 最適化する, "geometric": 最適化なし
        self.wpts = wpts  # 最適化の点数
        # 初期化
        self.set_mode(self.mode)
        self.reset_state()

    def _design_filters(self) -> None:
        f1, f2 = self.band
        if self.design == "fit":
            self.a0, self.a1 = _fit_allpass_pair(self.fs, f1, f2, self.nsec, self.wpts)
        elif self.design == "geometric":
            f0 = np.geomspace(f1, f2, self.nsec)
            f1s = np.clip(f0 * np.sqrt(2.0), f1 * 1.05, f2 / 1.05)
            self.a0 = np.array([_a_from_fc(fc, self.fs) for fc in f0], dtype=float)
            self.a1 = np.array([_a_from_fc(fc, self.fs) for fc in f1s], dtype=float)
        else:
            raise ValueError("design must be 'fit' or 'geometric'")
        # 係数（b,a）を生成
        self._ba0 = [(np.array([a, 1.0], float), np.array([1.0, a], float)) for a in self.a0]
        self._ba1 = [(np.array([a, 1.0], float), np.array([1.0, a], float)) for a in self.a1]

    def set_mode(self, mode: str) -> None:
        mode = mode.lower()
        if mode not in ("pos", "neg"):
            raise ValueError("mode must be 'pos' or 'neg'")
        self._sgn = +1.0 if mode == "pos" else -1.0
        self.mode = mode
        self._design_filters()

    def reset_state(self) -> None:
        self._zi0 = [np.zeros(1, dtype=complex) for _ in range(self.nsec)]
        self._zi1 = [np.zeros(1, dtype=complex) for _ in range(self.nsec)]

    def _cascade_lfilter(self, x: NDArray, ba_list, zi_list):
        y = x
        for i, (b, a) in enumerate(ba_list):
            y, zi_list[i] = signal.lfilter(b, a, y, zi=zi_list[i])
        return y

    def process(self, x: NDArray) -> NDArray:
        """フィルタ処理"""
        x = np.asarray(x).astype(complex, copy=False)
        y0 = self._cascade_lfilter(x, self._ba0, self._zi0)
        y1 = self._cascade_lfilter(x, self._ba1, self._zi1)
        return 0.5 * (y0 + (1j * self._sgn) * y1)

 実施例具体例 で示した複素信号 x(t) = 0.8 \ e^{1.5 \times 2\pi it} + 0.5 \ e^{6.0 \times 2\pi it} + 0.5 \ e^{-3.0 \times 2\pi it} に対して作成したフィルタをかけます。
# 信号作成
fs = 100.0
total_time = 10.
t = np.arange(int(total_time*fs)) / fs
x = 0.8*np.exp(1j*2*np.pi*1.5*t) + 0.5*np.exp(1j*2*np.pi*6.0*t) + 0.5*np.exp(-1j*2*np.pi*3.0*t)
x += 0.01*(np.random.default_rng(0).standard_normal(len(t)) + 1j*np.random.default_rng(1).standard_normal(len(t)))

# フィルタをかける (正の周波数)
flt = OneSidedIIR(fs, band=(0.3, 10.0), nsec=6, mode="pos", design="fit")
y_pos = flt.process(x)

# 負の周波数を抽出する
flt.set_mode("neg")
y_neg = flt.process(x)
周波数ごとの強度を求める関数は以下のように定義しました。
def spectrum(x, fs, nfft=8192):
    X = np.fft.fftshift(np.fft.fft(x, nfft))
    f = np.fft.fftshift(np.fft.fftfreq(nfft, d=1/fs))
    mag = np.abs(X) + 1e-12
    return f, 20*np.log10(mag/mag.max())
元の信号, 正の周波数のみになるようにフィルタをかけたもの, 負の周波数になるようにしたものは以下です。複素信号の実部を青色, 虚部をオレンジ色で表示しています。

正の周波数を抽出したものは実部が虚部に先行しており、e^{i\omega t} で \omega > 0 のように振る舞っていることが分かります。負の周波数を抽出したものは虚部が実部に先行しています。




周波数解析の結果は以下です (再掲)。元信号と正の周波数のフィルタ通過後の結果です。負の周波数が大きく減衰していることが分かります。

 参考文献以下を参考にしました。
基礎からの周波数分析（29）－「ヒルベルト変換と解析信号」
信号処理に欠かせないオールパスフィルタ―の概要・使い方・利用例
An Infinite Impulse Response (IIR) Hilbert Transformer Filter Design Technique for Audio