これは数理最適化 Advent Calendar 2022 の 23 日目の記事です。

はじめに

ディジタルフィルタとは、離散時間信号に対し所望の周波数帯域の成分だけを通過させて他は阻止する機能のことを言う。中でもIIR (Infinite Impulse Response) フィルタはフィードバック機構を用いることにより無限長の畳み込み計算を行うフィルタのことを言う。 IIRフィルタは特定の制約を満たしつつ所望の周波数特性が得られるよう適切に設計する必要がある。フィルタ設計方法として様々な方法が提案されており、その中には最適化手法を利用したものも数多くある。

本記事は主に参考文献[1]の1章の内容を基に、最適化によるIIRフィルタ設計を実装してその感触を掴んでみたという内容になっている。文献ではPSO（Particle Swarm Optimization）で最適化問題を解いているが、今回は手始めに（個人的に）実装ハードルが低い焼きなまし法による設計を試している。

ディジタルフィルタの基本

まずはディジタルフィルタについての基本事項を述べる。ここでは最低限のみ触れるので詳しくは参考文献[2][3]等を参照のこと。

フィルタの設計問題

IIRフィルタの処理は、離散時間の入力信号を $x [ i ]$ 、出力信号を $y [ i ] \ \ (i = 0, \cdots, \infty)$ とすると、以下のような差分方程式で表現できる。ただし、 $i \lt 0$ の場合は $x[i] =0, y[i ] = 0$ とする。

$\displaystyle y [i ] = -\sum_{m = 1}^{M} b_m y [i - m ] + \sum_{n = 0}^{N} a_n x [i - n ] \tag{1}$

この $b_m \in \mathbb{R}, (m = 1, \cdots, M)$ と $a_n \in \mathbb{R}, (n = 0, \cdots, N)$ がフィルタ係数で、所望の周波数特性が得られるようにこの係数を決めるのがフィルタの設計にあたる。

式（1）の右辺第1項は過去の出力値を利用するフィードバック機構になっており、第2項が過去の入力値を利用するフィードフォワード機構になっている。このフィードバック機構があることによってフィルタの設計が不適切な場合、出力信号が無限に発散してしまう可能性がある。このことを不安定といい、逆に発散する心配のないフィルタは安定という*1。従って安定かつ所望の周波数特性が得られるようにIIRフィルタを設計する必要がある。

Z変換

連続信号を周期 $T$ でサンプリングした離散信号はデルタ関数を用いて以下のように表すことができる。

$\displaystyle x_s(t) = \sum_{n = 0}^{\infty} x(nT) \delta(t - nT) \tag{2}$

このサンプリングされた信号をラプラス変換すると、

$\displaystyle \begin{align} \mathcal{L}[x_s(t)] &= \int_{0}^{\infty} \sum_{n = 0}^{\infty} x(nT) \delta(t - nT) \exp(-st) dt \\ &= \sum_{n = 0}^{\infty} x(nT) \exp(-nsT) \end{align} \tag{3}$

となる。ここで $z = \exp(sT)$ としたのがZ変換であり、次のような式になる*2。

$\displaystyle X(z) = \sum_{i = 0}^{\infty} x [ i ] z^{-i} \tag{4}$

ラプラス変換では $s = j\omega$ を代入することでフーリエ変換を得ていたが、Z変換も同様に $z = e^{j\omega}$ を代入すれば離散時間フーリエ変換を得ることができる（ $\omega$ は角周波数で周波数を $f$ とすると $\omega=2\pi f$ である）。そのため、この操作で周波数領域における特性を知ることができる。ただし、離散時間フーリエ変換ができるのは $\sum_{i = 0}^{\infty} | x [ i ]| \lt \infty$ という条件を満たす場合のみであるため、これを満たさない信号に対して $z = e^{j\omega}$ を代入する操作を行っても意味はない。

伝達関数

フィルタの周波数特性を見るには伝達関数が便利である。式（1）の両辺をZ変換にすると、

$\displaystyle Y (z) = -\sum_{m = 1}^{M} b_m Y(z) z^{-m} + \sum_{n = 0}^{N} a_n X (z) z^{-n}\tag{5}$

となる。ここで $H(z) = Y(z)/X(z)$ は伝達関数と呼ばれ、式（5）を整理してIIRフィルタの伝達関数を求めると以下のようになる。

$\displaystyle H (z) =\frac{Y(z)}{X(z)} = \frac{\sum_{n = 0}^{N} a_n z^{-n}}{1 + \sum_{m = 1}^{M} b_m z^{-m}} \tag{6}$

この伝達関数に対して、 $z = e^{j \omega}$ を代入すると伝達関数を離散時間フーリエ変換した結果が得られる。これを以下のように $H(\omega)$ と定義する。

$\displaystyle H (\omega) = \frac{\sum_{n = 0}^{N} a_n e^{-jn \omega}}{1 + \sum_{m = 1}^{M} b_m e^{-jm\omega}} \tag{7}$

$H(\omega)$ はフィルタの周波数特性であり、このフィルタによって各周波数成分がどれくらい減衰するかという振幅特性 $| H(\omega) |$ と、どれくらい位相が遅れるかという位相特性 $\angle H(\omega)$ がわかる*3。振幅特性はどの周波数成分を通過させてどの周波数成分を阻止するかを左右するのでフィルタにおいて当然重要な特性だが、位相特性も重要で通過させたい帯域に関してはなるべく線形であることが望ましい。位相特性が線形でないと、出力信号の波形が歪んでしまうためである。

極と安定性

伝達関数の分母=0を満たす点を極という。フィルタが安定であるためには全ての極の絶対値が $1$ 未満である必要がある。そのため、分母のフィルタ係数 $b_m$ は $1 + \sum_{m = 1}^{M} b_m z^{-m} = 0$ の解が全て絶対値 $1$ 未満となるように定めなければならない。このことは最適化問題における制約条件となる。

焼きなまし法の適用

定式化

まず所望の周波数特性を $H_D(\omega)$ とする。例えば、以下のような周波数特性を設定する。

$\displaystyle H_D(\omega) = \begin{cases} e^{-j \omega \tau_d} & \omega \in \Omega_p \\ 0 & \omega \in \Omega_c \end{cases} \tag{8}$

ここで、 $\tau_d$ は所望の群遅延、 $\Omega_p$ は通過させたい角周波数の集合、 $\Omega_c$ は阻止したい角周波数の集合である。この理想的な周波数特性に近づくよう式（6）、（7）の係数を決めたいので、係数を並べたベクトル $\mathbf{x} = [a_0, \cdots, a_N, b_1, \cdots, b_M ]$ を変数として、以下のような最適化問題を解くことになる。

$\displaystyle \begin{align} \min_{\mathbf{x}}& \ \ J(\mathbf{x})\\ \rm{s.t.}& \ \ | p_i(\mathbf{x}) | < 1 \ \ (i = 1, \cdots, M) \end{align} \tag{9}$

ここで $J(\cdot)$ は理想周波数特性との誤差関数、 $p_i(\mathbf{x})$ はフィルタ係数が $\mathbf{x}$ の時の極である。誤差関数は色々考えられるがよく使われるのは以下の2つ*4。

最大絶対誤差

$\displaystyle J(\mathbf{x}) =\max_{\omega \in \Omega} |H_D(\omega) - H(\omega; \mathbf{x})| \tag{10}$

平均2乗誤差

$\displaystyle J(\mathbf{x}) = \frac{1}{|\Omega|}\sum_{\omega \in \Omega} |H_D(\omega) - H(\omega; \mathbf{x})|^2 \tag{11}$

$\Omega$ は誤差計算に使用する周波数点の集合である。

さて、式（9）には極に関する制約条件があるがこのままでは扱いにくいので、制約条件違反をペナルティ項として組み込んだ新たな目的関数 $F(\mathbf{x})$ を作り、こちらを焼きなましで解くことにする。

$\displaystyle \begin{align} \min_{\mathbf{x}}& \ \ F(\mathbf{x}) = \min_{\mathbf{x}} \ \ (J(\mathbf{x}) + c \phi(\mathbf{x})) \\ \phi(\mathbf{x}) &= \begin{cases} \max_{i} |p_i(\mathbf{x})|^2 & \max_i | p_i(\mathbf{x}) | \geq 1 \\ 0 & \max_i | p_i(\mathbf{x}) | < 1 \end{cases} \end{align} \tag{12}$

$c$ はペナルティの強さを決めるハイパーパラメータになっている。

極の計算方法

式（12）のままだと更新の度に分母=0の $M$ 次方程式を解くなどして安定性を判断する必要がある。これでは効率が悪いため少し工夫する。まず簡単のためフィルタ係数の次数 $N, M$ は偶数とする*5。そうすると式（6）は次のように変形できる。

$\begin{align} \displaystyle H (z) &= \frac{\sum_{n = 0}^{N} a_n z^{-n}}{1 + \sum_{m = 1}^{M} b_m z^{-m}} \\ &= a0 \frac{\prod_{n = 1}^{N/2} (1 - q_n z^{-1})(1 - \bar{q}_n z^{-1})}{\prod_{m = 1}^{M/2} (1 - p_m z^{-1})(1 - \bar{p}_m z^{-1})} \\ &= a0 \frac{\prod_{n = 1}^{N/2} (1 - 2\rm{Re} (q_n) z^{-1} + |q_n|^{2} z^{-2}) } {\prod_{m = 1}^{M/2} (1 - 2\rm{Re}(p_m) z^{-1} + |p_m|^{2} z^{-2}) } \end{align} \tag{13}$

ここで、 $a_0 \in \mathbb{R}, p_m, \bar{p}_m, q_n, \bar{q}_n \in \mathbb{C}$ であり、 $\bar{p}_m, \bar{q}_n$ は $p_m, q_n$ の複素共役である。こう変形するとシステムの極は $p_m, \bar{p}_m$ なので直接的にわかる。従って、 $\mathbf{x} = [a_0, p_1, \cdots, p_{M/2}, q_1, \cdots, q_{N/2} ]$ *6のように直接極を変数として扱う形で式（12）の最適化問題を解けば更新毎の安定性の判別は簡単にできる。

最適化アルゴリズム

ここでは実際に実装したアルゴリズムの流れを述べる。焼きなまし法自体の解説は割愛するので、詳しくは参考文献[5]等を参照。

まず、理想周波数特性の関数 $H_D(\cdot)$ 、誤差計算に用いる周波数点集合 $\Omega$ 、フィルタ次数 $M, N$ 、ペナルティの強さ $c$ 、繰り返し回数 $L$ 、更新量の標準偏差 $\sigma$ 、初期温度 $T$ 、冷却率 $\alpha$ 、使用する誤差関数 $J(\cdot)$ はアルゴリズムのパラメータとして入力する。そして以下の処理を行う。

$\mathbf{x} = [a_0, p_1, \cdots, p_{M/2}, q_1, \cdots, q_{N/2} ]$ の各要素を正規分布 $\mathcal{N}(0, 1)$ に従う乱数で初期化する。
$\mathbf{x}' = \mathbf{x} + \Delta \mathbf{x}$ を求める。ここで $\Delta \mathbf{x}$ は各要素が正規分布 $\mathcal{N}(0, \sigma^2)$ に従う乱数。
$F(\mathbf{x}') \lt F(\mathbf{x} )$ であれば、 $\mathbf{x} \leftarrow \mathbf{x}'$ に更新する（ $F(\cdot)$ は式（12）の目的関数）。 $F(\mathbf{x}') \geq F(\mathbf{x})$ の場合は、 $\exp ( \frac{F(\mathbf{x}) - F(\mathbf{x}')}{T})$ の確率で $\mathbf{x} \leftarrow \mathbf{x}'$ に更新する。
温度を $T \leftarrow \alpha T$ として更新する。
繰り返し回数が $L$ 回未満なら2.に戻る。 $L$ 回に達したら、探索した中で最も $F(\mathbf{x})$ が小さかった $\mathbf{x}$ を出力して終了。

実装と実験

ここまで述べたことを実装して実験してみた。なお、ここに示す実験用コードはgithubにも上げている。

github.com

焼きなまし法の実装と実験

焼きなまし処理のコード

まずは焼きなまし法でフィルタ係数を求める処理の実装を以下に示す。

import pickle
import numpy as np
import matplotlib.pyplot as plt

def compute_mean_squared_error(H_, ref_h, omega_list):
    h = np.zeros(len(omega_list), dtype=np.complex)
    for i, omg in enumerate(omega_list):
        h[i] = H_(omg)

    e = np.real((ref_h - h) * np.conjugate(ref_h - h))
    return np.mean(e)

def compute_abs_max_error(H_, ref_h, omega_list):
    h = np.zeros(len(omega_list), dtype=np.complex)
    for i, omg in enumerate(omega_list):
        h[i] = H_(omg)
    e = np.abs(ref_h - h)
    return np.max(e)

# フィルタの周波数特性を計算する関数
def frequency_characteristic_func(a0, p_array, q_array, omega):
    denomi = 1.0
    for p in  p_array:
        denomi *= (1 - p * np.exp(-1j * omega)) *\
            (1 - np.conjugate(p) * np.exp(-1j * omega))

    numer = 1.0
    for z in q_array:
        numer *= (1 - z * np.exp(-1j * omega)) *\
            (1 - np.conjugate(z) * np.exp(-1j * omega))
    return a0 * (numer/denomi)

def obj_func(p_array, q_array, a0, ref_h, omega_list, c, error_func):
    H = lambda omega: frequency_characteristic_func(a0, p_array, q_array, omega)
    max_p2 = np.max(np.real(p_array * np.conjugate(p_array)))
    e = error_func(H, ref_h, omega_list)

    return e + (c * max_p2 if max_p2 >= 1.0 else 0.0)

# 焼きなまし法のメイン処理
def annealing(HD, omega_list, M, N, c, L, sigma, T, alpha, error_func,
              init_a0=None, init_p_array=None, init_q_array=None):

    # 理想周波数特性の関数から各周波数点に対する周波数特性を計算する
    h = np.zeros(len(omega_list), dtype=np.complex)
    for i, omg in enumerate(omega_list):
        h[i] = HD(omg)

    # 初期化
    if init_a0 is None:
        a0 = 1.0
    else:
        a0 = init_a0
    if init_p_array is None:
        p_array = np.random.randn(M//2) + 1j * np.random.randn(M//2)
    else:
        p_array = np.copy(init_p_array)
    if init_q_array is None:
        q_array = np.random.randn(N//2) + 1j * np.random.randn(N//2)
    else:
        q_array = np.copy(init_q_array)

    # 目的関数F
    F = lambda p_array_, q_array_, a0_: \
        obj_func(p_array_, q_array_, a0_, h, omega_list, c, error_func)

    now_cost = F(p_array, q_array, a0)

    # 最終的に出力するフィルタ係数を初期化
    best_p_array = np.copy(p_array)
    best_q_array = np.copy(q_array)
    best_a0 = a0
    best_cost = now_cost

    for i in range(L):
        # ロールバック用に保持
        save_a0 = a0
        save_p_array = np.copy(p_array)
        save_q_array = np.copy(q_array)

        a0 += np.random.randn() * sigma
        p_array += (np.random.randn(M//2) + 1j * np.random.randn(M//2)) * sigma
        q_array += (np.random.randn(N//2) + 1j * np.random.randn(N//2)) * sigma

        tmp_cost = F(p_array, q_array, a0)
        if np.exp((now_cost - tmp_cost)/T) > np.random.rand():
            now_cost = tmp_cost
            print("now_cost=", now_cost, "T=", T, "i=", i)
            if now_cost < best_cost:
                best_a0 = a0
                best_p_array = np.copy(p_array)
                best_q_array = np.copy(q_array)
                best_cost = now_cost
                print("best_cost = ", best_cost)
        else:
            # ロールバック
            a0 = save_a0
            p_array = np.copy(save_p_array)
            q_array = np.copy(save_q_array)
        T *= alpha

    return best_a0, best_p_array, best_q_array, best_cost

基本的には先述のアルゴリズムに書いた通りの実装になっている。annealingという関数が実際の焼きなましの処理を行う関数である。 $F(\cdot)$ の計算などを工夫すればもっと効率化できるかもしれない。

実験

それでは実際に理想的な周波数特性を設定し、この焼きなまし処理を用いてフィルタ係数を求めてみる。理想的な周波数特性は式（8）の形式において、通過域の集合を $\Omega_p = \{\omega |0.4 \pi \leq \omega \leq 0.6 \pi \}$ 、阻止域の集合を $\Omega_c = \{\omega | 0 \leq \omega \leq 0.3 \pi, 0.7 \pi \leq \omega \leq \pi \}$ 、群遅延を $\tau_d = 12$ と設定して実験した。なお、 $\Omega_p, \Omega_c$ を見てわかる通り、急激な変化を避けるため通過域と阻止域の間に特に周波数特性を指定しない遷移域も設けている。この理想周波数特性のグラフは以下のようになる。

焼きなましのパラメータについては、

$M = 8, N = 8, c = 5, L = 10000, \sigma = 0.01, T=10, \alpha = 0.998$

とした。誤差関数については両方の方法を試している。

また、問題の性質上、得られる結果が初期値や乱数に依存して大きく変わるので、この実装では初期値を変えながら100回最適化処理を行い、最も $F(\mathbf{x})$ が小さくなったものを採用している。

このような実験を行うコードは以下のようになっている。

def visualize_characteristic(H, omega_list, save_file_name):
    h = np.zeros(len(omega_list), dtype=np.complex)
    for i, omg in enumerate(omega_list):
        h[i] = H(omg)

    fig = plt.figure()
    amp_h = np.abs(h)
    ax1 = fig.add_subplot(2, 1, 1)
    ax1.plot(omega_list, amp_h, marker=".")
    ax1.set_ylim(-0.1, 1.1)
    ax1.set_ylabel("amplitude")

    angle_h = np.angle(h)
    ax2 = fig.add_subplot(2, 1, 2)
    ax2.plot(omega_list, angle_h, marker=".")
    ax2.set_ylim(-np.pi - 0.2, np.pi + 0.2)
    ax2.set_xlabel("omega")
    ax2.set_ylabel("phase")
    fig.savefig(save_file_name)
    plt.close()

if __name__ == '__main__':
    np.random.seed(0)

    # 理想的な周波数特性の関数
    HD = lambda omg: np.exp(- 1j * 12 * omg) \
        if omg >= 2 * 0.2 * np.pi and omg <= 2 * 0.3 * np.pi  else 0.0

    # 誤差関数に利用する誤差計算に使用する周波数点の集合
    omega_list1 = np.linspace(0, 2 * 0.15 * np.pi, 200)
    omega_list2 = np.linspace(2 * 0.2 * np.pi, 2 * 0.3 * np.pi, 200)
    omega_list3 = np.linspace(2 * 0.35 * np.pi, np.pi, 200)
    omega_list = np.concatenate([omega_list1, omega_list2, omega_list3])

   # 理想的な周波数特性のグラフを作成
    visualize_characteristic(HD, omega_list, "ideal.png")

    # 100回焼きなましをして最もよかったものを採用
    best_cost = 1e9
    for i in range(100):
        a0_, p_array_, q_array_, cost =\
            annealing(HD, omega_list, 8, 8, 5,
                      10000, 0.01, 10, 0.998,
                      compute_mean_squared_error) # 実験する誤差関数を指定
        if cost < best_cost:
            a0 = a0_
            p_array = np.copy(p_array_)
            q_array = np.copy(q_array_)
            best_cost = cost

    # 求めたフィルタの周波数特性の関数
    H = lambda omega: frequency_characteristic_func(a0, p_array, q_array, omega)

    # 求めたフィルタの周波数特性のグラフを作成
    visualize_characteristic(H, np.arange(0, np.pi, 0.01), "solution.png")

    # 制約条件を満たしていることを確認
    for p in p_array:
        print("|p|=", np.abs(p))

    print("final cost=", best_cost)

    # 求めたフィルタ係数を保存
    save_dict = {"a0":a0, "p_array":p_array, "q_array":q_array}
    with open('filter_coef.pickle', 'wb') as f:
        pickle.dump(save_dict, f)

これを実行して得たIIRフィルタの周波数特性は以下のようになった。

誤差関数が最大絶対誤差の場合

誤差関数が平均2乗誤差の場合

目視の印象だが通過域については最大絶対誤差で求めた方が平均2乗誤差に比べて理想の周波数特性に近いように見える。一方で、阻止域については最大絶対誤差だと振幅特性が0付近まで落ち切っておらず平均2乗誤差の方が優位に見える。

フィルタ処理の実装と実験

フィルタ処理のコード

得られたフィルタ係数を用いて実際の離散信号に対してフィルタをかけてみる。式（13）を展開して式（6）の形式にした上で式（1）の計算でフィルタをかけても良いが、ここでは式（13）の形式のまま2次のフィードバック（もしくはフィードフォワード）処理を縦続に接続することでフィルタ処理を実現している。つまり、2次フィルタの出力をそのまま次の2次フィルタに入力するような形にしている。これについては実装を見た方がわかりやすいだろうと思うので早速実装を貼る。

class IIRFilter:
    def __init__(self, a0_, p_array_, q_array_):
        self.a0 = a0_
        self.p_array = np.copy(p_array_)
        self.q_array = np.copy(q_array_)

    def _feedback_filter(self, x, p):
        b1 = -np.real(p) * 2
        b2 = np.real(p * np.conjugate(p))
        y = np.zeros(x.shape[0])
        for i in range(x.shape[0]):
            s = x[i]
            if i - 1 >= 0:
                s += -b1 * y[i - 1]
            if i - 2 >= 0:
                s += -b2 * y[i - 2]
            y[i] = s
        return y

    def _feedforward_filter(self, x, q):
        a1 = -np.real(q) * 2
        a2 = np.real(q * np.conjugate(q))
        y = np.zeros(x.shape[0])
        for i in range(x.shape[0]):
            s = x[i]
            if i - 1 >= 0:
                s += a1 * x[i - 1]
            if i - 2 >= 0:
                s += a2 * x[i - 2]
            y[i] = s
        return y

    def filtering(self, x):
        y = np.copy(x)
        for q in self.q_array:
            y = self._feedforward_filter(y, q)
        for p in self.p_array:
            y = self._feedback_filter(y, p)
        y *= self.a0
        return y

実験

この実装を使っていくつかの信号にフィルタをかけてみる。この実験用コードについては上述のgithubを参照。

$\cos (0.1 \pi t)$

まずは低周波数帯の遮断域にある周波数の信号をフィルタリングしてみる。以下のように減衰された信号が得られるが、誤差関数が最大絶対誤差の場合は平均2乗誤差の場合に比べて減衰が甘い。この結果は上述の振幅特性通りの結果と言える。

$\sin (0.4 \pi t)$

次に通過域ギリギリにある周波数の信号を試す。誤差関数が最大絶対誤差の場合は $t=0$ 付近以外ほとんど減衰のない出力が得られているが、平均2乗誤差では少し減衰した結果になっている。

$\sin (1.0 \pi t)$

高周波数帯の遮断域にある周波数の信号。低周波数帯の時と同様に、誤差関数が最大絶対誤差の場合は減衰しきれていない結果になっている。

$\cos (0.1 \pi t) + \sin (0.4 \pi t) + \sin (1.0 \pi t)$

最後に混合した信号で試す。通過域の成分のみ抽出できているように見える。波形としてきれいなのは、遮断域をきちんと遮断できている平均2乗誤差の方か。

所感

上でも述べたが良い解が得られるかどうかは初期値と乱数にかなり依存し、あまり好ましくない局所最適解が得られるケースも多々あった。特に誤差関数として最大絶対誤差を選択した場合では、振幅特性がほとんど0になるような解が得られるケースも多かった。そのため今回は初期値を変えながら何度も焼きなましを行う方法で対処したが、そもそも焼きなましで解くには非線形性が強すぎる問題だった可能性がある。そう考えると、PSOのように複数の個体を有する最適化手法を使うのは妥当な選択なのかもしれない。この辺は次の記事辺りでも比較してみたいところ。

また、フィルタの次数 $M, N$ を増やすほど最適化が難しくなる。変数が増えるので当然と言えば当然だが、今回の試行では $M=8, N=8$ （9変数）程度がまともな解が得られる限界だったのは少し残念だった。