ari23の研究ノート

メーカ勤務エンジニアの技術ブログです

MENU
トップ > Tech > ディジタル信号処理｜SOS型IIRフィルタ

ディジタル信号処理｜SOS型IIRフィルタ

Tech DSP Python

ScipyでIIRフィルタを設計すると、SOS型というものが出てきます。今回はこのSOS型IIRフィルタについて整理します🐜

IIRフィルタのおさらい

SOS型の説明の前に、一度IIRフィルタを復習します。

IIRフィルタとは、比較的低い次数で急峻な振幅特性が得られますが、フィルタを通すと発散してしまう不安定性をもつ可能性があることが特徴です。

ところで、IIRフィルタを差分方程式で書くと次の通りです。1

\displaystyle y[n] = \sum_{k=1}^{N}a_{k}y[n-k] + \sum_{k=0}^{M}b_{k}x[n-k]

a_{0} を導入し、式変形します。

\displaystyle \begin{align} a_{0}y[n] &= \sum_{k=1}^{N}a_{k}y[n-k] + \sum_{k=0}^{M}b_{k}x[n-k] \\ \sum_{k=0}^{M}b_{k}x[n-k] &= a_{0}y[n] - \sum_{k=1}^{N}a_{k}y[n-k] \end{align}

a の符号を入れ替えて置き直します。（ k=1 から k=0 となっていることに注意）

\displaystyle \sum_{k=0}^{M}b_{k}x[n-k] = \sum_{k=0}^{N}a_{k}y[n-k]

ここでZ変換を導入します。（ディジタル信号処理 ｜FIRとIIRの違いを参照）

\displaystyle \sum_{k=0}^{N} a_{k} z^{-k} Y(z) = \sum_{k=0}^{M} b_{k} z^{-k} X(z)

伝達関数 H(z) を導入すると、以下のように書けます。

\displaystyle H(z) = \frac{Y(z)}{X(z)} = \frac{\sum_{k=0}^{M} b_{k} z^{-k}} {\sum_{k=0}^{N} a_{k} z^{-k}}

この伝達関数の分母がゼロ（根）になるとき=極で安定性を見ているのが、ディジタル信号処理｜IIRフィルタ設計(scipy.signal.iirdesign)でした。

SOS型とは

IIRフィルタの特徴として、不安定性を有すると述べましたが、その傾向は次数が大きくなればなるほど顕著に現れます。 これを防ぐため、複数の2次(N=M=2)のIIRフィルタ（BiQuadフィルタ、双2次フィルタともいう）で構成し直したものを、SOS(second order sections)型と呼びます。

式とイメージ図を以下に示します。

\displaystyle \begin{align} H(z) &= \frac{\sum_{k=0}^{M} b_{k} z^{-k}} {\sum_{k=0}^{N} a_{k} z^{-k}} \\ &= \left( \frac{b_{10}+b_{11}z^{-1}+b_{12}z^{-2}}{a_{10}+a_{11}z^{-1}+a_{12}z^{-2}} \right) \left( \frac{b_{20}+b_{21}z^{-1}+b_{22}z^{-2}}{a_{20}+a_{21}z^{-1}+a_{22}z^{-2}} \right) ... \end{align}

SOS型IIRフィルタ
SOS型IIRフィルタ

高次元のIIRフィルタは、SOS型にすると安定性を有する傾向にあります。以降でこれを確認してみましょう。

SOSの安定性

ここでは、scipy.signal.sosfiltのドキュメントにある例を使って、SOSの安定性を確認します。

開発環境

サンプルプログラムを実行する開発環境は次の通りです。Python3であれば動くと思います。

項目 内容
OS Windows 10
Python 3.6.8
Numpy 1.18.5
Scipy 1.4.1
Matplotlib 3.2.1

サンプルプログラム

サンプルプログラム(sample_tf_sos.py)は以下の通りです。自由にコピペして頂いて構いません。ぜひ実際に動かしてみてください。

#!/usr/bin/python3
# -*- coding: utf-8 -*-
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

def ax_plot_poles(ax, poles):
    for p in poles:
        x, y = np.real(p), np.imag(p)
        if np.sqrt(x**2+y**2) > 1:
            ax.plot(x, y, marker='x', c='r')
        else:
            ax.plot(x, y, marker='x', c='k')
    return ax

def main():
    b, a = signal.ellip(13, 0.009, 80, 0.05, output='ba')
    sos = signal.ellip(13, 0.009, 80, 0.05, output='sos')

    theta = np.linspace(0, 2*np.pi, 360)
    x_unit_circle = np.cos(theta)
    y_unit_circle = np.sin(theta)

    x_impulse = signal.unit_impulse(1000)
    y_ba = signal.lfilter(b, a, x_impulse)
    y_sos = signal.sosfilt(sos, x_impulse)

    fig = plt.figure()

    ax = fig.add_subplot(2, 2, 1)
    ax.plot(x_unit_circle, y_unit_circle, c='k', ls=':')
    ax_plot_poles(ax, signal.tf2zpk(b, a)[1])
    ax.axis('equal')
    ax.set_ylabel('ba')

    ax = fig.add_subplot(2, 2, 2)
    # ax.plot(x_impulse, label='Unit impulse')
    ax.plot(y_ba, label='ba')

    ax = fig.add_subplot(2, 2, 3)
    ax.plot(x_unit_circle, y_unit_circle, c='k', ls=':')
    ax_plot_poles(ax, signal.sos2zpk(sos)[1])
    ax.axis('equal')
    ax.set_ylabel('sos')

    ax = fig.add_subplot(2, 2, 4)
    # ax.plot(x_impulse, label='Unit impulse')
    ax.plot(y_sos, label='sos')

    plt.show()

if __name__ == '__main__':
    main()

実行結果は以下の通りです。

極配置とインパルス応答の比較
極配置とインパルス応答の比較

解説

上記サンプルプログラムでは、次数13の楕円フィルタを設計し、その安定性を通常型（tfまたはba）と、SOS型で評価しています。

設計値はどちらも同じですが、SOS型にすると安定性を有することがわかります。

SOS型は安定性が高いと言えますが、一方で双2次フィルタを多段で組むため計算量が増えるデメリットもあるので、IIRフィルタを設計する際にはよく考える必要があります。

使用している関数などについては、ディジタル信号処理｜IIRフィルタ設計(scipy.signal.iirdesign)を参照してください。

おわりに

IIRフィルタのSOS型についてまとめました。

IIRフィルタを設計すると、安定性がなかなか担保されないこともあるので、このSOS型をうまく使いたいですね。

参考になれば幸いです(^^)

参考文献

  1. 備忘録の意味合いが強く、変数の説明は省略してます。↩