Main Content

IIR フィルターの設計

IIR フィルターと FIR フィルター

IIR フィルターが FIR フィルターに勝る点は、一般的に対応する FIT フィルターよりも少ないフィルター次数で与えられた仕様を満たせることにあります。IIR フィルターは非線形位相をもちますが、MATLAB® ソフトウェアでのデータ処理は一般にオフラインで行われ、データ シーケンス全体がフィルター処理の前に使用可能な状態になっています。これにより、関数 filtfilt を使用した非因果的なゼロ位相フィルター処理が可能になり、IIR フィルターの非線形位相ひずみを除去できます。

標準的 IIR フィルター

標準的な IIR フィルターであるバタワース、チェビシェフ I 型、チェビシェフ II 型、楕円、およびベッセルの各フィルターはすべて、異なる方法で理想的な「ブリック ウォール (急峻なロールオフ特性)」をもつフィルターを近似します。

このツールボックスには、アナログ領域とデジタル領域の両方で (アナログでのみサポートされるベッセルを除く)、また、ローパス、ハイパス、バンドパス、およびバンドストップの構成で、これらすべての標準的な IIR フィルターを作成する関数が用意されています。また、ほとんどのフィルター タイプに対して、通過帯域の減衰量や阻止帯域の減衰量、および遷移幅で与えられたフィルター仕様に適合する最小のフィルター次数を求めることもできます。

その他の IIR フィルター

直接型フィルター設計関数 yulewalk では、指定した周波数応答関数を近似する振幅応答をもつフィルターが求められます。これは、マルチバンドの帯域フィルターを作成する 1 つの方法です。

また、パラメトリック モデリング、または、システム同定関数を使用して、IIR フィルターを設計することもできます。これらの関数は、パラメトリック モデリングで説明します。

汎用バタワース設計関数 maxflat は、汎用バタワース フィルターの設計の節で説明します。

IIR フィルター方法の概要

次の表に、このツールボックスで用意されているさまざまなフィルター方法をまとめ、これらの方法の実装に使用できる関数を記載します。

ツールボックスのフィルター方法と使用可能な関数

フィルター方法説明フィルター関数

アナログ プロトタイピング

連続 (ラプラス) 領域で標準的なローパス プロトタイプ フィルターの極と零点を使用して、周波数変換とフィルター離散化によりデジタル フィルターを実装

完全設計関数: besselfbuttercheby1cheby2ellip

次数推定関数: buttordcheb1ordcheb2ordellipord

ローパス アナログ プロトタイプ関数: besselapbuttapcheb1apcheb2apellipap

周波数変換関数: lp2bplp2bslp2hplp2lp

フィルター離散化関数: bilinearimpinvar

直接設計

区分線形振幅応答を近似することで、離散時間領域において直接デジタル フィルターを設計

yulewalk

汎用バタワース設計

極よりも零点の数が多いローパス バタワース フィルターを設計

maxflat

パラメトリック モデリング

設定された時間領域応答、または周波数領域応答を近似するデジタル フィルターが求められます (パラメトリック モデリング ツールについては、System Identification Toolbox™ のドキュメンテーションを参照してください)。

時間領域モデリング関数: lpcpronystmcb

周波数領域モデリング関数: invfreqsinvfreqz

アナログ プロトタイプ作成による標準的な IIR フィルター設計

このツールボックスで提供されている主な IIR デジタル フィルターの設計法は、標準的なローパス アナログ フィルターを等価なデジタル フィルターに変換することをベースとしています。以下の節では、フィルターの設計法と、サポートされているフィルター タイプの特性の概要を説明します。フィルター設計プロセスの詳細な手順は、IIR フィルター設計での特別トピックを参照してください。

標準的な完全 IIR フィルターの設計

フィルター設計関数を使用して、ローパス、ハイパス、バンドパス、または、バンドストップ構成の、任意の次数のフィルターを簡単に作成できます。

フィルター設計関数

フィルター タイプ

設計関数

ベッセル (アナログのみ)

[b,a] = besself(n,Wn,options)

[z,p,k] = besself(n,Wn,options)

[A,B,C,D] = besself(n,Wn,options)

バタワース

[b,a] = butter(n,Wn,options)

[z,p,k] = butter(n,Wn,options)

[A,B,C,D] = butter(n,Wn,options)

チェビシェフ I 型

[b,a] = cheby1(n,Rp,Wn,options)

[z,p,k] = cheby1(n,Rp,Wn,options)

[A,B,C,D] = cheby1(n,Rp,Wn,options)

チェビシェフ II 型

[b,a] = cheby2(n,Rs,Wn,options)

[z,p,k] = cheby2(n,Rs,Wn,options)

[A,B,C,D] = cheby2(n,Rs,Wn,options)

楕円

[b,a] = ellip(n,Rp,Rs,Wn,options)

[z,p,k] = ellip(n,Rp,Rs,Wn,options)

[A,B,C,D] = ellip(n,Rp,Rs,Wn,options)

既定では、これらの関数からはそれぞれローパス フィルターが返されます。正規化周波数 (ナイキスト周波数 = 1 Hz) を使用して、希望するカットオフ周波数 Wn を設定するだけです。ハイパス フィルターでは、関数のパラメーターのリストに 'high' を追加します。バンドパスまたはバンドストップ フィルターの場合は、Wn を通過帯域のエッジ周波数を含む 2 要素ベクトルとして指定します。バンドストップ構成に 'stop' を追加します。

次に、デジタル フィルターの設計例をいくつか示します。

[b,a] = butter(5,0.4);                    % Lowpass Butterworth
[b,a] = cheby1(4,1,[0.4 0.7]);            % Bandpass Chebyshev Type I
[b,a] = cheby2(6,60,0.8,'high');          % Highpass Chebyshev Type II
[b,a] = ellip(3,1,60,[0.4 0.7],'stop');   % Bandstop elliptic

シミュレーション用などにアナログ フィルターを設計するには、最後に 's' を追加して、カットオフ周波数をラジアン/サンプル単位で指定します。

[b,a] = butter(5,0.4,'s');      % Analog Butterworth filter

すべてのフィルター設計関数では、出力引数の数に応じて、伝達関数、零点-極-ゲイン、または、状態空間線形システムのモデル表現でフィルターが返されます。一般的に、丸め誤差による数値的な問題が生じる可能性があるため、伝達関数の使用は避けてください。代わりに、zp2sos を使用して 2 次セクションに変換できる零点-極-ゲイン形式を使用し、その後に、変換した 2 次セクションを使用してフィルターを解析または実装します。

メモ

すべての標準的 IIR ローパス フィルターは、極端に低いカットオフ周波数には適していません。このため、非常に狭い通過帯域をもつローパス IIR フィルターを設計する代わりに、通過帯域の広いフィルターを設計して、入力信号を間引く方法を使用します。

周波数領域仕様に対する IIR フィルターの設計

このツールボックスには、与えられた必要条件のセットを満たす最小のフィルター次数を計算する次数選択関数が用意されています。

フィルター タイプ

次数推定関数

バタワース

[n,Wn] = buttord(Wp,Ws,Rp,Rs)

チェビシェフ I 型

[n,Wn] = cheb1ord(Wp,Ws,Rp,Rs)

チェビシェフ II 型

[n,Wn] = cheb2ord(Wp,Ws,Rp,Rs)

楕円

[n,Wn] = ellipord(Wp,Ws,Rp,Rs)

これらの関数は、フィルター設計関数と合わせて使うと便利です。たとえば、1,000 Hz から 2,000 Hz が通過帯域、その両側から 500 Hz 離れた位置に阻止帯域があり、サンプリング周波数が 10 kHz、通過帯域でのリップルが 1 dB 以内で、阻止帯域の減衰量が少なくとも 60 dB のバンドパス フィルターを設計するとします。関数 butter を使用してこれらの仕様を満たすには、以下のようにします。

[n,Wn] = buttord([1000 2000]/5000,[500 2500]/5000,1,60)
[b,a] = butter(n,Wn);
n =
    12
Wn =
    0.1951    0.4080

同じ必要条件を満たす楕円フィルターは、次のように設計します。

[n,Wn] = ellipord([1000 2000]/5000,[500 2500]/5000,1,60)
[b,a] = ellip(n,1,60,Wn);
n =
    5
Wn =
    0.2000    0.4000

これらの関数は、アナログ フィルターの場合に加え、ほかの標準帯域構成も扱うことができます。

標準的な IIR フィルター タイプの比較

このツールボックスには 5 つのタイプの標準的な IIR フィルターが用意されており、それぞれ異なる状況に適しています。本節では、各フィルター タイプの基本的なアナログ プロトタイプを示し、主な特性の概要を説明します。

バタワース フィルター

バタワース フィルターでは、アナログ周波数 Ω = 0 および Ω = ∞ において理想的なローパス フィルター応答に対する最善のテイラー級数近似が提供されます。任意の次数 N について、振幅二乗応答は、これらの位置で 2N – 1 個のゼロ微分をもちます (Ω = 0 および Ω = ∞ において "最大フラット")。応答は全体的に単調で、Ω = 0 から Ω = ∞ に滑らかに減少します。Ω = 1 においては |H(jΩ)|=1/2 です。

チェビシェフ I 型フィルター

チェビシェフ I 型フィルターでは、通過帯域内に Rp dB の等リップルを組み込むことにより、通過帯域全体にわたって理想の周波数応答と実際の周波数応答との差の絶対値が最小になるようにします。阻止帯域応答は、最大フラットとなります。通過帯域から阻止帯域への遷移はバタワース フィルターより急激です。Ω = 1 においては |H(jΩ)|=10Rp/20 です。

チェビシェフ II 型フィルター

チェビシェフ II 型フィルターでは、阻止帯域内に Rs dB の等リップルを組み込むことにより、阻止帯域全体にわたって理想の周波数応答と実際の周波数応答との差の絶対値が最小になるようにします。通過帯域応答は、最大フラットとなります。

阻止帯域は、I 型フィルターほど速くゼロに近づきません (また、偶数値のフィルター次数 n では、まったくゼロに近づきません)。ただし、多くの場合、通過帯域にリップルがないことは重要な利点になります。Ω = 1 においては |H(jΩ)|=10Rs/20 です。

楕円フィルター

楕円フィルターは通過帯域と阻止帯域の両方で等リップルとなります。一般的にこれは、サポートされる任意のフィルター タイプにおいて、その最小次数で必要条件を満たします。フィルターの次数 n、通過帯域リップル Rp デシベルおよび阻止帯域リップル Rs デシベルが与えられると、楕円フィルターによって遷移幅が最小化されます。Ω = 1 においては |H(jΩ)|=10Rp/20 です。

ベッセル フィルター

アナログ ベッセル ローパス フィルターでは、ゼロ周波数において群遅延が最大フラットになり、通過帯域全体にわたってほぼ一定の群遅延が維持されます。このためフィルター処理された信号は、通過帯域周波数範囲内で元の波形を維持します。アナログ ベッセル ローパス フィルターが周波数マッピングを介してデジタル フィルターに変換されると、この最大フラットという特性は失われます。Signal Processing Toolbox™ は、アナログの場合のみベッセル フィルター設計関数を完全にサポートします。

ベッセル フィルターは一般に、阻止帯域において十分な減衰を行うために、他のフィルターより高い次数を必要とします。Ω = 1 においては |H(jΩ)|<1/2 となり、フィルター次数 n が増加するにつれて減少します。

メモ

上記のローパス フィルターは、アナログ プロトタイプ関数 besselapbuttapcheb1apcheb2ap および ellipap を使用して作成されています。これらの関数により、カットオフ周波数 1 rad/s をもつ適切なタイプの n 次アナログ フィルターの零点、極、ゲインが求められます。完全なフィルター設計関数 (besselfbuttercheby1cheby2ellip) は、設計プロセスの最初の段階でプロトタイピング関数を呼び出します。詳細については、IIR フィルター設計での特別トピックを参照してください。

同様のプロットを作成するには、n = 5 とし、必要に応じて、Rp = 0.5 および Rs = 20 を使用します。たとえば、楕円フィルターのプロットを作成するには、次のようにします。

[z,p,k] = ellipap(5,0.5,20);
w = logspace(-1,1,1000);
h = freqs(k*poly(z),poly(p),w);
semilogx(w,abs(h)), grid
xlabel('Frequency (rad/s)')
ylabel('Magnitude')

直接形 IIR フィルターの設計

このツールボックスで "直接設計法" と言う場合は、離散領域での仕様に基づいてフィルターを設計する IIR 設計の手法を指しています。アナログ プロトタイプ作成法とは異なり、直接設計法は標準的なローパス、ハイパス、バンドパス、バンドストップ構成に制限されません。これらの関数では、任意の、場合によってはマルチバンドの周波数応答をもつフィルターが設計されます。本節では、特にフィルター設計を目的とした関数 yulewalk について説明します。パラメトリック モデリングでは、Prony 法、線形予測法、スティグリッツ・マクブライド法、逆周波数設計など、直接型ともみなされ得るその他の手法について説明します。

関数 yulewalk では、指定した周波数応答に適合するように再帰型 IIR デジタル フィルターが設計されます。yulewalk という関数名は、フィルターの分母係数を求める方法から付けられたものです。つまり、この関数では指定した理想的な振幅二乗応答の逆 FFT が求められ、その結果として得られた自己相関関数のサンプルを使用して、拡張ユール・ウォーカー方程式を解けます。ステートメント

[b,a] = yulewalk(n,f,m)

周波数の振幅特性がベクトル fm で与えられたものに近似する、n 次の IIR フィルターの分子と分母の係数 n+1 個を含む行ベクトル ba が返されます。f は 0 ~ 1 の範囲の周波数点を表すベクトルで、ここで、1 はナイキスト周波数を表します。m は、f 点において指定した振幅応答を含むベクトルです。fm は、マルチバンド応答を含む、任意の区分線形形状の振幅応答を表すことができます。この関数に相当する FIR 関数は fir2 で、この関数も任意の区分線形振幅応答に基づいたフィルターを設計します。詳細については、FIR フィルターの設計を参照してください。

yulewalk には位相情報を入力することができません。また、結果として得られるフィルターの最適性に関しては、保証をするものではありません。

yulewalk を使用してマルチバンド フィルターを設計し、指定する周波数応答と実際の周波数応答をプロットします。

m = [0   0   1   1   0   0   1   1   0 0];
f = [0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 1];
[b,a] = yulewalk(10,f,m);
[h,w] = freqz(b,a,128)
plot(f,m,w/pi,abs(h))

汎用バタワース フィルターの設計

ツールボックスの関数 maxflat を使用すると、汎用バタワース フィルター、すなわち、零点と極の数が異なるバタワース フィルターを設計できます。これは、零点より極の方が計算時間がかかる実装では望ましいといえます。maxflat は、1 つの次数のみでなく、分子と分母に対して 1 つずつ、2 つの次数を設定できるという点を除けば、関数 butter と同じです。これらのフィルターは "最大フラット" となります。これは、結果として得られるフィルターが任意の分子および分母の次数に対して最適であり、0 およびナイキスト周波数 ω = π における導関数の最大数がともに 0 であることを意味します。

たとえば、2 つの次数が等しいとき、maxflatbutter と同じです。

[b,a] = maxflat(3,3,0.25)
b =
    0.0317    0.0951    0.0951    0.0317
a =
    1.0000   -1.4590    0.9104   -0.1978
[b,a] = butter(3,0.25)
b =
    0.0317    0.0951    0.0951    0.0317
a =
    1.0000   -1.4590    0.9104   -0.1978

ただし、maxflat では極より多くの零点をもつフィルターを設計できるため、用途が広いといえます。

[b,a] = maxflat(3,1,0.25)
b =
    0.0950    0.2849    0.2849    0.0950
a =
    1.0000   -0.2402

maxflat に対する 3 番目の入力は、周波数 0 ~ 1 で 1/2 の目標振幅応答をもつ、"電力半値周波数" です。

また、'sym' オプションを使用して、最大フラットの特性をもつ線形位相フィルターを設計することもできます。

maxflat(4,'sym',0.3)
ans =
    0.0331    0.2500    0.4337    0.2500    0.0331

maxflat アルゴリズムの詳細については、Selesnick と Burrus による [2] の参考文献を参照してください。