ケィオスの時系列解析メモランダム

時系列解析、生体情報解析などをやわらかく語ります

【デジタルフィルタの基礎5】sinc関数で作るローパスFIRフィルタ:理想と現実のはざま

教科書や文献には「理想のローパスフィルタ」というものが登場します。もちろん「理想」は人それぞれ違ってかまわないのですが、通常ここで言う理想とは、「ある周波数(カットオフ周波数 f _ c)より低い成分は完全に通し、それより高い成分は完全に遮断する」という性質をもつフィルタを指します。このフィルタの周波数応答をグラフで描くと、下図左のような、急に立ち上がった「壁」のような形になります。このような特性を、「理想矩形(くけい)特性」や「ブリックウォール(レンガの壁)フィルタ」と呼びます。

理想のローパスフィルタ(左)とそのインパルス応答(右)

 周波数領域での形が決まれば、それをフーリエ変換することで時間領域のインパルス応答を計算できます。つまり、フィルタの係数列を求めることができます。したがって、「理想のローパスフィルタ」が周波数領域で矩形関数ならば、そのインパルス応答は上図右のような sinc 関数になります。
 とはいえ、ここでは、FIRフィルタとしてローパス特性を実現したいと考えています。したがって、ここでは、sinc 関数が時間方向に無限の幅を持つ(無限サポートである)ことが問題になります。FIRフィルタとして実装するためには、インパルス応答を有限の幅(有限サポート)にしなければなりません。
 そこで考えられる方法が、無限に広がる sinc 関数を「どこかで切って」有限長にするという変形です。最も単純な方法は、ある幅だけを残してそれ以外をゼロにする、いわゆる「切り取り(トランケーション)」です。しかし、この方法では時間領域で急に切り落とすため、周波数領域ではギブズ現象と呼ばれる「リップル(ripple: 波打ち)」が生じ、理想とは違った特性が現れてしまいます (下図)。
ローパスフィルタの構造。フィルタの構造を表す用語については、この図を参照してください。

 そこで登場するのが窓関数(window function)です。窓関数は、sinc 関数を単にバサッと切るのではなく、両端をなめらかにゼロへ近づけるように重みをかけることで、周波数領域のリップルを抑えつつ有限長化を行うための道具です。ハミング窓・ハニング窓・ブラックマン窓など、さまざまな窓が知られており、それぞれが、「リップルの抑制の程度」や「遷移帯域の鋭さ」といった性質のバランスを持っています。
 このように、周波数領域では矩形特性(=理想)を思い描きつつ、時間領域では sinc を有限長化するための現実的な工夫が必要になる——これが、FIR ローパスフィルタ設計における「理想と現実のはざま」です。窓関数法はその橋渡しをするもっとも基本的な手法であり、多くのデジタル信号処理の教科書で最初に紹介されています。
 今回は、sinc関数を下絵として、ローパスFIRフィルタを作る基礎的知識について解説します。

1. 理想のローパスフィルタ:周波数領域では矩形、時間領域ではsinc

 離散時間信号 \displaystyle{
x[n]
} に対して、理想的なローパスフィルタの周波数応答 \displaystyle{
H _ {\mathrm{d}}(f)
} は、

\displaystyle{
H_{\mathrm{d}}(f)=
\begin{cases}
1 & (|f|\le f_c) \\
0 & (|f|>f_c)
\end{cases}
}

のような 完全な「長方形(矩形)」 になります(\displaystyle{
f _ c
} はカットオフ周波数)。

理想的なローパスフィルタ

  • パスバンド(\displaystyle{|f|\le f _ c})では 完全に 1(=そのまま通す)
  • ストップバンド(\displaystyle{|f|\gt f _ c})では 完全に 0(=完全カット)
  • 遷移帯域はゼロ(急激にストンと落ちる)

という、教科書で登場する “理想的なフィルタ” です。この「矩形関数」を逆離散時間フーリエ変換(逆DTFT) すると、

\displaystyle{
h_{\mathrm{d}}[n]
= \int_{-1/2}^{1/2} H_{\mathrm{d}}(f)\, e^{j 2\pi f n}\, df
= \int_{-f_c}^{f_c} e^{j 2\pi f n}\, df 
}

となり、

\displaystyle{
h_{\mathrm{d}}[n] =
\begin{cases}
\dfrac{\sin(2\pi f_c n)}{\pi n} & (n\neq 0) \\
2 f_c & (n = 0)
\end{cases}
}

がえられます。これは、sinc 型のインパルス応答です。 標準的な sinc 関数

\displaystyle{
\mathrm{sinc}(x)=\frac{\sin(\pi x)}{\pi x}
}

を用いると、上の式は

\displaystyle{
h_{\mathrm{d}}[n]
= 2 f_c \, \mathrm{sinc}(2 f_c n)
}

のように書くことができます。

\displaystyle{f _ c = 0.05} のときの、\displaystyle{ h _ {\mathrm{d}}[n] = 2 f _ c \, \mathrm{sinc}(2 f _ c n) } のプロット。

2. 理想 sinc から現実的なフィルタをデザイン

 理想インパルス応答 \displaystyle{
h _ {\mathrm{d}}[n]
} には、実装上、次の大きな問題があります。それは、
 無限長である
ということです。係数が無限に続くと計算が面倒なので、ここから、「sinc をうまく切り取って有限長FIRにする」ということを考えます。

2-1 sinc関数の両端をぶった切って捨てる

最も素朴なアイデアは、

中心付近だけ \displaystyle{2M+1} 点だけ残して、あとは 0 にしてしまう

というものです。
 つまり、長さ 2M+1矩形窓(rectangular window)

\displaystyle{
w_{\text{rect}}[n-M] =
\begin{cases}
1 & (n = 0,1,\dots,2M)\\
0 & \text{それ以外}
\end{cases}
}

を使って、

\displaystyle{
h[n] = h_{\mathrm{d}}[n-M]\, w_{\text{rect}}[n-M], \quad n = 0,1,\dots,2M
}

のように、切り取ってやれば、有限長になります。私は、生体信号のオフライン処理を想定しているので、非因果であることは気にしません (因果にしたければシフトしてください)。
 このとき、フィルタ係数は具体的には

\displaystyle{
h[n] =
\begin{cases}
\dfrac{\sin\bigl(\omega_c (n-M)\bigr)}{\pi (n-M)} & (n \neq M)\\
\dfrac{\omega_c}{\pi} & (n = M)
\end{cases}
}

です(\displaystyle{
n=0,\dots,2M
})。
 ぶった切るというのは簡単ですが、ローパスフィルタを作りたいので、周波数応答はどうなるか?が本質的な関心事です。
 下の図は、フィルタ長 N(あるいは半長 M)を変化させたときの 周波数応答(左)インパルス応答(右) を示しています。右図の青線がインパルス応答で、薄い灰色は sinc関数そのものです (インパルス応答の合計を1にしているので両者の高さは異なります)。左図の青実線がインパルス応答に対応する周波数応答です。

sinc関数を有限幅でトランケートしたしたときの周波数応答。トランケートする幅を変化させた場合。

 上の図左の周波数応答には、N をある程度大きくしても、次のような Gibbs 現象 が現れます。

  • カットオフ付近に、ギブス現象特有の「つの状のオーバーシュート(でっぱり)」が生じる
  • パスバンドおよびストップバンドに リンギング(リップル) が発生する
  • 遷移帯域が理想よりも広がる

これは、理想的な sinc 形インパルス応答を矩形窓で有限長に切り出すと必ず現れる、典型的な挙動です。
 これでも「ある程度は満足できる」という人もいるでしょうが、もう少し良い特性を得たい。すなわち

  • オーバーシュート(つの状のはみ出し)を抑えたい
  • パスバンド/ストップバンドのリップルをもっと小さくしたい
  • 遷移帯域をもう少し鋭くしたい

と考える人も多いはずです。
 そこで次に登場するのが、窓関数(window function) を用いた改良手法です。

2-2 インパルス応答の端をするっと0に近づける:窓関数法

 矩形窓で突然バッサリと sinc を切り落とすのではなく、なだらかに減衰する窓関数を掛けてトランケートすることで、周波数応答の急激な不連続を和らげ、Gibbs 現象を大幅に抑制することができます。
 典型的には、以下のような窓を用いることで特性が改善します。

  • Hann(ハン)窓
  • Hamming(ハミング)窓
  • Blackman(ブラックマン)窓
  • Kaiser(カイザー)窓

 これらの窓関数は、時間領域で滑らかにゼロへ向かう形をしており、その分だけ周波数領域での振る舞いも滑らかになります。結果として、

  • リップルが減少し、
  • オーバーシュートが小さくなり、
  • ストップバンド減衰が改善する

というメリットが得られます。
 ここでは,最もよく使われる窓の一つである Hamming(ハミング)窓 を用いた例を見てみます。Hamming 窓は、矩形窓のように端で突然ゼロになるのではなく、端に向かってゆるやかに減衰する形をしており、その結果として周波数応答のリンギング(リップル)を大幅に低減できます。
 Hamming 窓  w[n] は、フィルタ長を  N とすると、次式で定義されます。

\displaystyle{
w[n] = 0.54 - 0.46 \cos\!\left( \frac{2\pi n}{N-1} \right),
\qquad n = 0,1,\dots, N-1.
}

この窓を理想インパルス応答 \displaystyle{
h _ d[n]
} に掛け合わせ、

\displaystyle{
h[n] = h_d[n] \, w[n]
}

とすることで、有限長の FIR フィルタが得られます。
 下の図は、「矩形窓」を次第に「Hamming 窓」に変形した場合に、周波数応答(左)とインパルス応答(右)がどのようになるかを描いたものです。右図の桃色破線が、窓関数の形を示しています。ただし、高さが実際とは異なることに注意して下さい (実際は最大値が1になります)。

矩形窓からハミング窓へ変形すると周波数応答の振動が弱まる

上図左の周波数応答を見ればわかるように、窓関数を両端でスムーズに落ち込む Hamming 窓 に変えれば、

  • パスバンドおよびストップバンドのリップルが明らかに減少し、

  • カットオフ付近のオーバーシュート(ギブス現象)が大きく抑えられ、

  • ストップバンド減衰が向上する

といった改善が視覚的に確認できます。
 さらに、トランジションバンドの幅を狭くして、理想の周波数応答に近づけたいのであれば、下の図のようにフィルタ長を長くする必要があります。このフィルタ長については、分析したいデータの長さとのせめぎ合いです。

窓の幅を広げれば理想のフィルタに近づく

3. まとめ:まずは“イメージ”をつかめ

 今回は、理想ローパスフィルタ(矩形特性)を出発点にして、sinc を有限長に切り出して FIR フィルタを作るという、一連の流れをできるだけ直感的につかめるよう図を使って説明しました。
 実際のフィルタ設計では、窓関数を選ぶ・フィルタ長を決める・遷移帯域やストップバンド減衰の仕様を満たす――といった技術的な調整が必要になりますが、今回は詳細には立ち入りません。
 実際の応用では、Hamming 窓以外にも Hann、Blackman、Kaiser、Chebyshev など多くの窓関数があり、 目的に応じて「リップルをどこまで許すか」「遷移帯域をどれだけ鋭くしたいか」などを考えながら選択する必要があります。
 ここで得た直感を出発点として、ぜひ読者自身の目的に応じて学習を深め、設計の自由度を広げていってください。

※ もし記事の中で「ここ違うよ」という点や気になるところがあれば、気軽に指摘していただけると助かります。質問や「このテーマも取り上げてほしい」といったリクエストも大歓迎です。必ず対応するとは約束できませんが、できるだけ今後の記事で扱いたいと思います。それと、下のはてなブログランキングはあまり信用できる指標ではなさそうですが (私のブログを読んでいる人は、実際とても少ないです)、押してもらえるとシンプルに励みになります。気が向いたときにポチッとしていただけたら嬉しいです。