直感的に理解する信号処理とフィルタ【STEP2：インパルス応答とフィルタ】

# はじめに

・キッチンスケールで小麦粉の重さを量る
・Bluetoothイヤホンで音楽を聴く
・GPSで今いる場所を確認する

これら全ての基礎になる技術が**信号処理**です。
そして信号処理で扱うあらゆる信号には影のようにつきまとってくるものがあります。

そう、**ノイズ**ですね。 

　・電源ラインから入り込んでくるノイズ
　・モータの振動によるノイズ
　・電子回路のクロストークノイズ
 
原因は様々ですが通常、信号とノイズは必ずセットで扱わなければなりません。

最も有効なノイズ対策はノイズ源そのものを除去することですが、例えば工場での重量計測において、計測中は周囲で動いてる装置を全部止めて振動を無くしましょう！なんて現実的ではありません。

そこで活躍するのがフィルタ技術、とりわけデジタル(ディジタル)フィルタです。

このシリーズではフィルタ設計法習得に向け、フィルタ設計と信号処理の基礎に関する直感的な理解を目標に解説していきたいと思います。

:::note warn
なるべく数式を排して取っつきやすい説明を行っていくつもりなので、やや不自然な言い回しや誤解を招きかねない表現があるかもしれません。
そういった箇所を見つけられた場合はすみませんが、ご指摘頂けると幸いです。
:::

本記事はその **STEP2：インパルス応答とフィルタ** です。

他STEPへのリンクは下記です↓
[STEP:0 導入編](https://qiita.com/otter1602/items/826ee58a45dda366f9b2)
[STEP:1 フィルタ特性の見方](https://qiita.com/otter1602/items/20ae3138a8ef3530df7d)
[STEP:3 周波数領域における伝達関数](https://qiita.com/otter1602/items/a9292b9d91c00dc7273b)
[STEP:4 ノイズ除去実践編](https://qiita.com/otter1602/items/9719c0850f6cffc2d006)

簡単なフィルタ設計&性能評価をやってみたい方は下記もどうぞ↓
[FIRローパスフィルタを作ってみよう！](https://qiita.com/otter1602/items/83ec804e822f0842d701)

# 1. 前回のおさらい
前回の記事で、デジタルフィルタの周波数特性について解説しました。

周波数特性は **振幅特性**と **位相特性** に分けられ、それぞれフィルタへ入力した信号の大きさと(単位円上で見た時の)位置をどう変化させて出力するか、を表すものでした。

・振幅特性
![Step1_FIR_Amp.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4103962/6e56c2d4-5697-4b8c-8fb3-2ec0f5b36c6f.png)

・位相特性
![Step1_FIR_Phase.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4103962/84f66c61-d642-4311-a9c3-0c8243efac38.png)

・フィルタからの出力信号：振幅は0になり、位相は+90degされる
![Step1_+90deg_24.375Hz.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4103962/59d423de-32d4-40ec-a652-6ae19c940fcb.png)

さて、位相特性について前回は **線形位相特性** というものを解説しました。

・線形位相特性(先程の位相特性をアンラップしたもの)
![Step1_FIR_Phase_Unwrap.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4103962/fea1a794-b128-45d8-8fc8-96802ece14de.png)

このフィルタの場合は20Hzまでの信号を入力すると、入力信号の周波数と位相の変化が線形の関係になるのでした。

では、線形位相特性があるならば、 **非線形位相特性** もあるのでしょうか？
答えはYESです。

# 2. 非線形位相特性
いきなりですが、下記にフィルタの周波数特性を示します。

・振幅特性
![Step2_IIR_Amp.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4103962/a7cd4c12-1c66-47f4-a744-4c0bc6eea2c3.png)

・位相特性
![Step2_IIR_Phase.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4103962/18bbddf0-69a7-43f2-bd92-c5d2d2036f03.png)

こちらのフィルタですが前回示したフィルタと同様、20Hz以上の信号を入力するとゲインが0になる振幅特性を持っています。

今回注目したいのは位相特性の方です。上記に示したのはラップされた位相特性ですが、明らかに前回示したフィルタの位相特性とは異なっています。

アンラップした位相特性を確認してみましょう。

・位相特性(アンラップ)
![Step2_IIR_Phase_Unwrap.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4103962/6081f007-0245-4f09-8ae4-e791e65b6578.png)

・再掲：前回示したフィルタの位相特性(アンラップ)

![Step1_FIR_Phase_Unwrap.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4103962/b7cd20cf-3c2e-47e7-9b02-c412655ec8c0.png)

アンラップした位相特性を見ると分かりやすいかと思いますが、今回のフィルタは線形位相特性を持っていません。

すなわち、**非線形位相特性** を持つフィルタになります。

前回説明した通り、線形位相特性を持つフィルタでは入力信号の周波数に対する出力信号の位相変化の関係が完全に直線になります。

しかし、非線形位相特性を持つフィルタではその関係は成り立たず、位相変化は周波数ごとにバラバラです。そのため、フィルタを通した結果として **「波形が歪む」** という現象が起こります。(「ゆがむ」ではなく「ひずむ」です！)

例を出しましょう。

こちらに非線形位相特性を持つフィルタに矩形波を入力し、出力信号を取得したグラフを示します。青線が入力信号の矩形波、赤線が出力信号です。

![Step2_iir_phase_distortion_static.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4103962/727b82de-564a-42c9-9559-e424616d7508.png)

ご覧の通り、出力信号では入力信号には無かった**ツノ**が立っていますね。

ご存知と思いますが、矩形波は基本波とその高調波の足し合わせでできています。
ということは、元の矩形波は様々な周波数を含んだ信号です。

なので非線形位相特性を持つフィルタを通した結果、矩形波に含まれる周波数ごとに位相がバラバラに変化して出力され、位相変化の大きい高調波が **ツノ**になって表れています。 

比較のため、線形位相特性を持つフィルタで同様のことをやってみます。
今度は青線が入力信号の矩形波、緑線が出力信号です。

![Step2_fir_phase_nodistortion_static.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4103962/e74b96fd-476b-4af2-84cf-461aaa7c7d17.png)

どうでしょうか？ **ツノ**は無く、高調波の一部が除去されて正弦波の形になっていますね。 
すなわち、出力信号の波形は元の波形の形を保ったまま、振幅特性に従った振幅のみカットされたということになります。

さて、それではフィルタが線形位相特性を持つか、それとも非線形位相特性を持つか、一体どこで決まるのでしょうか？

最も簡単なのは**FIRフィルタならば線形位相特性、IIRフィルタならば非線形位相特性** という分け方です。

# 3. FIRフィルタ/IIRフィルタとは
デジタル(ディジタル)フィルタには、大きく分けて下記の2種類があります。

(1) FIRフィルタ：**Finite Impulse Response Filter = 有限インパルス応答フィルタ** 
(2) IIRフィルタ： **Infinite Impulse Response Filter = 無限インパルス応答フィルタ** 

さっそくそれぞれの説明に入りたいところですが、まずインパルス信号とインパルス応答について簡単に解説します。

## ・インパルス信号
下記図を見てください。図のように時刻 n=0のときだけ"1"、その他は0になるような信号を**インパルス信号** と呼びます。

![Step2_Impulse_Signal.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4103962/e2479fa0-4dab-4c3b-9c49-b656c7d43232.png)


 このインパルス信号を何かしらの **「システム」** に入力したとき、そのシステムがどんな出力信号を出すか、言い換えるとシステムがどんな応答を返すかを **インパルス応答** と呼びます。


※私の感覚ですが、「システム」という呼び方は制御工学とか信号処理の呼び方ですね。
　"入力信号->システム->出力信号" この関係はちょうど "引数->関数(メソッド)->戻り値"と
 　 同じか少なくとも似た関係です。

 インパルス応答がなぜ重要かというと、理想的にはインパルス信号には直流成分から超高周波な成分まで、あらゆる周波数成分が等しい大きさで含まれているからです。
 
 入力信号にあらゆる周波数成分が含まれているので、当然出力信号もあらゆる周波数成分に対する出力になりますね。
 つまり、インパルス応答が分かる = 各周波数における出力特性が分かるということです。

 そして、ここまで見てきたフィルタも「システム」の一種なので、下記が成り立ちますね。
 **「フィルタのインパルス応答 = 各周波数における出力特性 = フィルタの周波数特性」** 
 ※フーリエ変換/逆変換の関係のため、実際はイコールではありません。
 　あくまで、それぞれの持つ情報という視点では同じと言っても良いだろう、程度です。

 さて、そんな大事なインパルス応答ですが、インパルス信号を入力したシステムのふるまいによって、応答が継続する時間が変化します。

継続する時間が有限であれば、**有限インパルス応答** 、無限であれば**無限インパルス応答** と呼ばれます。といってもこれだけでは何が何やらなので、それぞれ図示します。

## ・有限インパルス応答

![Step2_Impulse_FIR.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4103962/8e0cf4b1-6c15-4e05-b2f3-734176911841.png)


この場合はインパルス応答はn=3で0に収束していますね。

## ・無限インパルス応答

![Step2_Impulse_IIR.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4103962/68ec5188-82b7-4a65-b4fe-bceb4d60973e.png)


この場合インパルス応答は小さくなっていくものの0には収束せず、n=20でもまだ続いています。実用上、計算機の扱える桁数には制限があるためいずれは0として扱われることになりますが、理論上は無限に続きます。

具体的なイメージのためひとつ例え話をしましょう。お寺の鐘をつくことを考えてみます。
一回だけ鐘をつくのがインパルス信号で、鐘の音の響きがインパルス応答です。

<有限インパルス応答>
一回鐘をついて、ボ～ンと音が響きます。でも、放っておくといずれ音は止みますね。
これがまさに有限インパルス応答です。

![Step2_FIR_Ex.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4103962/8f66ba8a-dcdc-4c53-95bc-449ae8d4628e.png)



<無限インパルス応答>
一回鐘をつくのは変わりません。鐘をついて、ボ～ンと音が響きます。
ですが今度はあらかじめ鐘の中にスピーカーか何かが置いてあって、直前の鐘の音を自動的に再生し続けるような仕掛けがされているような状況を考えます。
すると一回鐘をついた後は再生した音が勝手に鐘を響かせて、鐘をついていなくても音が
段々小さくなりつつ響き続ける、という状態になります。
これが無限インパルス応答です。

![Step2_IIR_Ex.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4103962/54ad4968-5219-40a9-9260-b8562269b145.png)

簡単に、と言いつつ長くなってしまいましたが、以上がインパルス応答に関する説明です。実はこれでFIRフィルタとIIRフィルタについてもほぼ説明したことになります。

なぜなら、フィルタの設計をする = フィルタの周波数特性を決めることですが、フィルタの周波数特性を決めることは、すなわち望むインパルス応答をするシステムを設計する、ということに他ならないからです。

なので本節の最初の疑問である、FIRフィルタとは何か？IIRフィルタとは？に対する答えは、次の通りです。

**(1) FIRフィルタ：インパルス信号を入力したとき、応答が有限で終了するフィルタ** 
**(2) IIRフィルタ：インパルス信号を入力したとき、応答が無限に継続するフィルタ** 

前節でFIRフィルタは線形位相特性、IIRフィルタは非線形位相特性を持つと述べました。
位相特性とインパルス応答の関係を見るためには、もう少し説明が必要になります。

# 4. 伝達関数
前節でインパルス応答について見ていきましたが、線形/非線形位相特性のようなフィルタの特性とインパルス応答を結びつけるために必要な概念である**伝達関数** について解説します。

## ・入力の平均化

少しだけ数式を使います。

たった今取得した入力信号の値を$z^0$、1回前に取得した値を$z^{-1}$と置き、さらに2回前に取得した値を$z^{-2}$と置きます。

これを一般化すると、$n$回前に取得した値は$z^{-n}$と置けます。
$z^{-n}$を使うと、$n=N$回前までの値の総和は下記のようになりますね。

$$
N回前までの値の総和 = \sum_{n=0}^{N} z^{-n}
$$

さて、突然ですが**移動平均** というものを皆さんご存知だと思います。ちょっとデータにノイズがのってるかな？とか、なんとなく波形を滑らかにしたいかな？なんて時にとても便利ですよね。
では、ちょっと考えて頂きたいんですが、先程示した総和の式を用いてN回移動平均の式はどう書けるでしょうか？
↓
↓
↓
↓
↓
↓
↓
↓
↓
↓
↓
↓
$$
N回移動平均 = \frac{1}{N}\sum_{n=0}^{N} z^{-n}
$$

答えはもちろん、求めた$n=N$回前までの値の総和を$N$で割ったものになります。
さらにこの$\frac{1}{N}$は最後にまとめて$N$で割るのではなく、値取得時に$\frac{1}{N}$をかけ算しても同じことなので、N回移動平均は下記のように書くこともできます。
$$
N回移動平均 = \sum_{n=0}^{N} \frac{1}{N} z^{-n}
$$
ここでもし、N回取得する値のうち最新の値が一番信頼できて、値が古くなるにしたがってあまり信頼できないことが分かっているとしましょう。

その場合は全ての値に毎回$\frac{1}{N}$をかけるよりも、何回前の値かによってかけ算する値を変える方が結果は正確になりそうですよね。

つまり、1回前の値はかなり信頼できるので$\frac{5}{N}$、N回前の最も古い値は全然信頼できないので$\frac{0.2}{N}$といった具合で重み付けをするということです。

分かりやすくするため、元々$\frac{1}{N}$固定だった値をnごとに変化する値として、$b_n$と置いてみましょう。すると、重み付きのN回移動平均の式は下記のように書けますね。

$$
重み付きN回移動平均 = \sum_{n=0}^{N} b_n z^{-n}
$$

※もちろん$b_n=\frac{1}{N}$であれば、重み無しの通常のN回移動平均と全く同じです。

## ・出力のフィードバック
もう一歩踏み込んでみましょう。さきほど重み付きN回移動平均について考えてみましたが、入力信号にとんでもない外れ値があった場合を想像してみましょう。

いくら重み付き移動平均とはいえ、しばらくはその外れ値の影響は残りそうですよね。
ですが、早く外れ値の影響を取り除きたい！となった場合、どうすれば良いでしょうか？

重み付きN回移動平均の場合、入力信号のみを使用して出力を計算していました。
では、外れ値があった場合にそれを検知して出力の変化を抑える手段として、自身の過去の出力を使うことを考えてみてはどうでしょうか？

入力信号と同様の考え方ですが、今度は$n$回前に出力した値を$z^{-n}$と置いてみます。
$z^{-n}$を使うと、$n=M$回前までの値の総和は下記のようになりますね。
※注意点として、出力側なので$n=0$ではなく、$n=1$から始まっています。
$$
M回前までの値の総和(出力側) = \sum_{n=1}^{M} z^{-n}
$$

さて、出力信号から外れ値を検知するといっても、入力信号が少々大きかったからといって、あんまり敏感に反応してしまうと出力が不安定になりそうですよね。

なので、こちらも移動平均を使いましょう。入力側と全く同様です。
単純移動平均ではなく、重み付き移動平均で考えてみると、下記のような式ができますね。

※さきほどは変化する値として$b_n$を置きましたが、今回は$a_n$としています。
$$
重み付きM回移動平均(出力側) = \sum_{n=1}^{M} a_n z^{-n}
$$

ではこの出力側の重み付き移動平均を使って、入力側の外れ値を検知し、出力の急な変化を抑えるような式を考えてみましょう。
といっても、難しいことはありません。単純に引き算するだけです。

すなわち、下記のようになります。
$$
出力信号(フィードバック付き) = \sum_{n=0}^{N} b_n z^{-n} - \sum_{n=1}^{M} a_n z^{-n}
$$

せっかく入力信号と出力信号の関係式を2種類導き出したので、それぞれを入出力をきちんと含んだ形で書いておきましょうか。

ここまで「n回前の入出力」というのを$z^{-n}$で表現してきたので、それにならって入力信号を$X(z)$、出力信号を$Y(z)$で表現します。

・入力信号からのみ出力を計算する場合
　元の式は下記です。
 $$
重み付きN回移動平均 = \sum_{n=0}^{N} b_n z^{-n}
$$
　これを$X(z)$、$Y(z)$を使って表現すると、
$$
Y(z) = (\sum_{n=0}^{N} b_n z^{-n}) X(z)
$$
と書けますね。


一方で、
・出力のフィードバックと入力信号から出力を計算する場合
　元の式はすぐ上で書きましたが、再掲します。

 
 $$
出力信号(フィードバック付き) = \sum_{n=0}^{N} b_n z^{-n} - \sum_{n=1}^{M} a_n z^{-n}
$$
　これを$X(z)$、$Y(z)$を使って表現すると、
$$
Y(z) = (\sum_{n=0}^{N} b_n z^{-n})X(z) - (\sum_{n=1}^{M} a_n z^{-n})Y(z)
$$
となります。

## ・フィルタの入出力関係と伝達関数

すでにお察しかもしれませんが、ここで導いた2つの式がそれぞれFIRフィルタ、IIRフィルタにおける入出力の関係式になります。

**入力信号のみで出力が構成されるのがFIRフィルタ、入力信号に出力フィードバックを加えて構成されるのがIIRフィルタ** です。

さらに、前節でフィルタは「システム」と呼ぶことができ、「システム」はプログラムでいうところの関数と同じか、少なくとも似た働きをする、と述べました。

フィルタがプログラムでいう関数と似ているならば、フィルタの振る舞いは引数を受け取って戻り値を返すまでの間にある処理と考えることができそうです。

もっと言うとフィルタの場合は処理といっても論理演算ではなく全て数値演算なので、処理の中身は入力値と出力値の比で表すことができるはずですね。

そこで、入力値と出力値の比を$H(z)$と表すことにしましょう。さきほど入力信号を$X(z)$、出力信号を$Y(z)$で表現することにしたので、すなわち$H(z)$は下記のように表されます。
$$
H(z) = \frac{Y(z)}{X(z)}
$$

こういう入出力の比で表される$H(z)$を制御工学や信号処理では**伝達関数**と呼びます。 

それでは、FIRフィルタとIIRフィルタの伝達関数をそれぞれ求めてみましょう。

・FIRフィルタ
FIRフィルタの入出力関係はこうでした。
$$
Y(z) = (\sum_{n=0}^{N} b_n z^{-n}) X(z)
$$

伝達関数はそのまま式変形するだけで求められます。
$$
H(z) =  \frac{Y(z)}{X(z)} = \sum_{n=0}^{N} b_n z^{-n}
$$
結局、FIRフィルタの伝達関数は
$$
重み付きN回移動平均 = \sum_{n=0}^{N} b_n z^{-n} = H(z)
$$
重み付きN回移動平均と同じになりますね。
すなわち、FIRフィルタは入力信号のみを使って出力を決定するということが伝達関数からも分かります。


・IIRフィルタ
一方、IIRフィルタの入出力関係はこうです。
$$
Y(z) = (\sum_{n=0}^{N} b_n z^{-n})X(z) - (\sum_{n=1}^{M} a_n z^{-n})Y(z)
$$

これも伝達関数はそのまま式変形で求められます。一応順を追って変形してみます。
$$
Y(z) + (\sum_{n=1}^{M} a_n z^{-n})Y(z) = (\sum_{n=0}^{N} b_n z^{-n})X(z)
$$
移項して、$Y(z)$でくくります。
$$
Y(z)(1 + \sum_{n=1}^{M} a_n z^{-n}) = (\sum_{n=0}^{N} b_n z^{-n})X(z)
$$
さらに変形して、
$$
H(z) = \frac{Y(z)}{X(z)} = \frac{(\sum_{n=0}^{N} b_n z^{-n})}{(1 + \sum_{n=1}^{M} a_n z^{-n})} 
$$
これが、IIRフィルタの伝達関数になります。

FIRフィルタの伝達関数と比較すると一目瞭然ですが、過去の出力をフィードバックするのがIIRフィルタということが伝達関数から確認できます。

# まとめ
ここまでフィルタの種類とインパルス応答、加えて伝達関数について述べてきました。

インパルス応答はインパルス信号を入力信号とした時の出力のことでした。
そして、さきほど求めた伝達関数はインパルス信号に限らず、より一般化されたフィルタの入出力関係を表すものでした。

なので結局、インパルス応答と伝達関数は本質的には入出力関係という、同じものを指していると言えます。

さて、フィルタの伝達関数が分かったことで入力信号に対する出力信号を求めることができるようになりました。

ですが、今回説明したインパルス応答と伝達関数はどちらも時間軸上でのお話です。
一方、これまでフィルタの振幅・位相特性はどちらも周波数軸上で見てきたはずです。
なので、下記のような疑問を持たれるかもしれません。

・インパルス応答や伝達関数を周波数軸上で当てはめるにはどうすれば良いのか？
・前回から引き継いだ疑問として、フィルタの応答遅延はどうやって求められるのか？

今回は道具立ての意味合いが大きい記事で、どちらかというと理論に寄った内容でした。
次回の記事では上記疑問に対する答えとして、もう少し実践的な内容を解説していきます。