【信号処理・制御工学の基礎】線形差分方程式の表現まで

はじめに

以下の記事の続き的な記事です。需要の無さは理解していますが自分の勉強と復讐を兼ねて記事を書いていきます。

s0sem0y.hatenablog.com

前回の復習

話の骨子は、 $n=k$ のときのみ $δ[n-k]=1$ を取るインパルス信号を導入することで、任意の入力信号を

$x[n] = \sum_k x[k]δ[n-k]$

と表現できることから始まり（前回記事では、 $a[k]$ って別の文字を置いてましたが、実際のところ上記で良い）、線形性と時不変性から

$\mathcal L(x[n]) = \mathcal L \left( \sum_k x[k]δ[n-k] \right) =\sum_k x[k] \mathcal L(δ[n-k]) =\sum_k x[k] h[n-k]$

と変形したのでした（詳しくは前回記事を）。

結局前回の話はというと、線形時不変システム $\mathcal L$ のインパルス応答 $h[n]=\mathcal L(δ[n])$ を用いると、任意の入力信号 $x[n]$ に対するシステム $\mathcal L$ の応答 $y[n]$ が

$\displaystyle y[n] = \mathcal L(x[n]) = \sum_{k=-∞}^{∞} x[n-k]h[k]=\sum_{k=-∞}^{∞} x[k]h[n-k]$

で求められるという結論を見たということです。これは $x$ と $h$ については可換であることも最後にちょろっと見ました。

入力信号をシステムに入れることを考えた時には、私達はシステムがどのように信号に対して加工をするのかを把握しなければなりません。

しかし、線形時不変システムの場合は入力に対してどのように加工を行うのかの情報が、インパルス応答に全て含まれており、具体的に行われる加工は、入力信号に対してインパルス応答を畳み込むということで達成されるのです。

インパルス応答でのシステムの考え方

前回の記事で、システムとは一体なんでしょうか、関数みたいなものかもしれないしそうではないかもしれない、ということを述べました。線形時不変システムの場合は、システムの方もあたかも何らかの信号（インパルス応答が本質）だと思ってしまって良いのです。

インパルス応答を作る＝線形時不変システムを作る

とりあえず私達は、線形時不変システム $\mathcal L$ なるものを設計することと、インパルス応答 $h[n]$ を設計することを同一視することができるようになりました。

それはつまり、線形時不変システムが変形が必ず $h[n]$ で達成されると分かっているのであれば、「なんかよくわからんシステムを作ってから、それがどのようなインパルス応答を持つのかを調べ、そして実際に畳み込み演算を行う」なんてことをするのではなく、最初から入力信号に畳み込みたい信号を作れば良いのです。その信号を設計することは線形時不変システムを設計していることにほかなりません。

どうせ、

$\displaystyle y[n] = \sum_{k=-∞}^{∞} x[n-k]h[k]$

によって計算されるのですから、 $h[k]$ という信号の方を作ってしまえばいいんですね。言い換えてしまえば、入力信号になにか適当な信号を畳み込んだら、線形時不変システムによる変換を行っているのと同じなのです！（それがどんな効力を持っているのかは、結果を見てみないとわからないが！）

仮に以下のケースで実験してみます。入力は

$x[n] = \sin(2π(0.03 n)) + 0.2 \times \mathcal N(0,1)$

f:id:s0sem0y:20171024073706p:plain

コイツに対して、

$h[n] = \sum_{m=0}^4 0.2δ([n-m])$

を畳み込みます。小難しい書き方していますが、要するにindexが $0〜4$ で値が $0.2$ となり、他は全部 $0$ となっているような信号を考えています。

f:id:s0sem0y:20171024074302p:plain

ただ信号 $h[n]$ を畳み込んでみただけなのに、随分と第二項のノイズ成分が減ったように見えます（これは、いわゆる移動平均フィルタに相当するものを使ったことになる）。

もしも、自分がやりたい処理に合わせて $h[n]$ を上手く設計できたら便利ですよね。ただし、時間領域での設計というのは線形システムに対する詳細をあまり知ることが出来ません（それでも試行錯誤的に設計して、良いものができるなら全然OK）。このために様々なフィルタ設計法（例えば周波数領域での設計）が提案されているというわけです。

補足

時間を意識するならば上記のような書き方になります（が、今はとりあえず畳みこめば線形時不変システムを使っていることになる例を見たかったのでそうしなかった）。先程は横軸がサンプリング点数になっていたが、こちらでは横軸を時間[sec]で扱っています。ただ、このサンプリング周波数周りは基本的に人間が管理しなければなりません。

特にサンプリング定理は重要で、「１秒間に100点をサンプリングするならば、その離散信号が復元できるのは50Hzの信号まで」です。つまり、上記の実験に置いて、ノイズ項が50Hz以上の周期性を持っていたとしても、それを知るすべはありません（むしろ、折り返し雑音によって波形は汚される）。

またnumpyの畳み込みはあくまで演算をしてくれるだけで、これが時系列信号であるのかどうかなんて知ったこっちゃないので、mode='same'でやっていいのかは状況によると思います（因果性を考える場合、mode='full'で計算し、はみ出た右のindexを捨てるべきと思われます）。

リアルタイムで動くシステム

ここで非常に重要な概念である「因果性」を導入します。簡単に言えば、これはシステムがリアルタイムで動くか否かを決定しうるものです。

因果性を有するというのは、そのシステムのインパルス応答 $h[n]$ が $n \geq 0$ のみに値を持つ場合を指します。システム $\mathcal L$ にインパルス信号 $δ[n]$ を入力した場合の応答がインパルス応答でした。ところでこのインパルス応答は $n=0$ にのみ値を取りますから、システムは $n=0$ 時点でのみ、たった１回の入力を受け取るのです。

もしもシステムのインパルス応答が $n = 0$ で値を持っているのであれば、それはシステムが現在の入力（ $n=0$ でのインパルス）に依存することを意味しており、 $n ＞ 0$ で値を持っていれば過去の値に依存する（ $n＞0$ 以降も、 $n=0$ のインパルスの影響が残っている）ことを意味します。

では、インパルス応答が $n＜0$ に値を持つというのはどういうことでしょうか。入力は $n=0$ にしか行われていないのです。つまり、 $n＜0$ 地点で、その地点から見れば未来の事象である $n=0$ の影響を受け取っているようなシステムになってしまうのです。

通常、ディジタル信号処理においては因果性があるもののみを重点的に扱います。すなわち、システムは現在かそれよりも過去の情報のみを使って出力をするように設計されるのです。過去の情報は通常メモリに保存しておいて後で使い回すことになります。

一方で非因果的なシステムは未来の情報を必要とするために、例えば $y[n]$ の出力を決定するために、次の入力 $x[n+1]$ を待つ必要があります。システムが値を出力してくれるのは、その情報が必要だった時刻よりも幾分か後になるということです。

ただし、これは「非因果的なシステムが使えないものである」ということを意味しません。出力が幾分か遅れても良い場合は勿論のこと、例えば信号の解析をじっくり腰を据えてやりたい場合は非因果的なシステムでも構わないです。ただ、通常「制御」や「信号処理」の場では応答性というのは安定性と並び非常に重要な指標となるためシビアであるケースが多いのは事実です。

線形差分方程式によるフィルタ設計

なるほど線形時不変システムのことを考える場合は、因果性に注意しつつ、インパルス応答（事実上単に信号 $h[n]$ )を上手く作ってやれば良さそうです。

もしもそれが出来るならばいいでしょう、ところで、 $h[n]$ の長さはどれくらいにしましょうか。値はどのようにしましょうか。どうすればどのような出力が得られるのか、いちいち実験するしか方法が無いのでしょうか。

実はそんなことはありません。ちゃんと数理的に解析する方法があります。周波数解析への道を開きましょう（教科書だと最初からこっちでスタートすると思われますが）。

FIRフィルタとIIRフィルタ

先程まで述べてきたように、線形時不変システムというのはインパルス応答で特徴づけることができます。そうした場合には、インパルス応答が有限の時間で止まるタイプの物と、無限の時間まで続くタイプのものに分かれます。これらを表すための言葉として、それぞれ

・FIRフィルタ(finite impulse response filter)

・IIRフィルタ(infinite impulse response filter)

があります。私達はインパルス応答を上手く作れば線形時不変システムを作れたことになることを知っていますが、世の中にはIIRフィルタなるものも存在するのです。つまり無限に長く続くインパルス応答を持つ線形時不変システムがあるということです。

インパルス応答を設計して、畳み込みを計算したいところですが、インパルス応答が無限に続くのであれば計算が出来ません。ということはIIRフィルタというものは実用性が無い理論上の存在なのでしょうか。

これを解決するのが差分方程式によるフィルタの表現です（実は線形時不変システムというのは線形差分方程式としての側面もあるのです。というより、普通こっちがスタート。連続なら線形微分方程式）。

以下は、とあるIIRフィルタの線形差分方程式による表現です。これは因果的なシステムを考えているものとして、 $y[n] = 0　if　n＜0$ としておきます。

$y[n] = y[n-1] + x[n]$

さてこのようなシステムのインパルス応答はどんなふうになるでしょうか。インパルス応答は入力に $δ[n]$ (n=0でのみ1となる)というインパルス信号を使うのでした。

$y[0] = y[-1] + δ[0] =0+ 1 = 1$

$y[1] = y[0] + δ[1] = 1 + 0 = 1$

$y[2] = y[1] + δ[2] = 1 + 0 = 1$

...というわけで、一生 $1$ を出力し続けるようです。したがって、上記の差分方程式によって表現された線形時不変システムを、インパルス応答の畳み込みを計算するという方法によって実装することは事実上不可能であるということです。

ただ、私達はインパルス応答が線形時不変システムを表現してくれることを知っているだけで使わなければならないわけでもありません。ですからケースバイケースで、インパルス応答を設計するのか、線形差分方程式を設計するのかを変えていくということになります。

実際

$y[n] = y[n-1] + x[n]$

なるシステムを使いたければ、現在の入力と前回の出力を足し合わせれば良いだけの話で、これだけで既にシステムの設計は済んでいるのです。わざわざインパルス応答を求めてから畳み込もうとしなくてもいいのです(インパルス応答を使わないにしても、無限の畳み込みを許容すれば、ちゃんとインパルス応答の畳込みによって線形時不変システムが表現できる事実は揺るがない)。

次に、以下のシステムを考えてみましょう。実はこれは、前回の記事の最初の方に出現させたFIRフィルタです。もう一度インパルス応答を見てみましょう。

$y[n] = x[n-1] + x[n]$

これは先ほどとは違い、前回の入力と現在の入力の和を取るというものです。これのインパルス応答は

$y[0] = δ[-1] + δ[0] = 0 + 1 = 1$

$y[1] = δ[0] + δ[1] = 1 + 0 = 1$

$y[2] = δ[1] + δ[2] = 0 + 0 = 0$

となって、以降、出力は得られません。結局のところ、有限の過去までの入力に依存するシステムは、入力が無くなったらいずれは出力も無くなります。インパルス応答とは、 $n=0$ にのみ $1$ という値を加えた際の応答なわけで、入力にしか依存しないシステムは必ずFIRフィルタとなります。

一方で、自分自身の出力をもう一度入力にするタイプのフィルタはIIRフィルタとなります。

以下に一般的なフィルタの形式を書いておきます。

FIRフィルタ

$\displaystyle y[n] = \sum_{k=0}^{K-1} b[k] x[n-k]$

結局有限の長さ(ここでは $K$ )の $b[k]$ との畳み込みで計算されるものをFIRフィルタというわけです（FIRフィルタはそもそもインパルス応答が有限である。明らかにFIRフィルタとは、有限のインパルス応答が $h[n] = b[n\$ となっている線形時不変システムに相当する）。

IIRフィルタ

$\displaystyle y[n] = \frac{1}{a[0]} \left( -\sum_{j=1}^{J-1} a[j]y[n-j] +\sum_{k=0}^{K-1} b[k]x[n-k] \right)$

IIRフィルタではインパルス応答との畳み込みを考えると無限遠まで計算が続くことになります。したがって明示的にインパルス応答との畳み込みという形では現実的な計算は与えられません。これを、出力が過去数点分フィードバックされていると考えることで、有限の区間で表現できるようになる場合もあります。ですから現実的にはこれが出来る範囲でシステムを設計するということになります（これは、決して全ての線形時不変システムを漏れ無く網羅できるわけではない）。

ちなみに、本当は $\frac{1}{a[0]}$ の項はあってもなくても良い（他の係数に吸収してもらえる）のですが、次の話に進むために入れておきます。お気づきだと思われますが、 $y[n-j]$ の係数を全て $0$ だと思ってしまえば、FIRフィルタにもなります。

したがってこちらは、現実的に表現できる線形時不変システムの線形差分方程式による表現と言ってもいいでしょう。

線形差分方程式によるフィルタ処理

まず以下の表現は、

$\displaystyle y[n] = \frac{1}{a[0]} \left( -\sum_{j=1}^{J-1} a[j]y[n-j] +\sum_{k=0}^{K-1} b[k]x[n-k] \right)$

$\displaystyle a[0]y[n] = -\sum_{j=1}^{J-1} a[j]y[n-j] +\sum_{k=0}^{K-1} b[k]x[n-k]$

$\displaystyle \sum_{j=0}^{J-1} a[j]y[n-j] =\sum_{k=0}^{K-1} b[k]x[n-k]$

と変形できることも知っておきましょう。こうすると幾分か見た目が綺麗になりますね。通常はこのような形式で書かれることが多いです。とにもかくにも、係数 $a,b$ の列を作ってやればフィルタを設計してやったことになるわけです。

これはだいたいscipy.signalにかなりたくさん実装があるので、すぐにできます（numpyのconvolveを使う人は多分いないと思う）。上記の差分方程式において、

scipy.signal.lfilter

最初の例と同じように、ノイズの乗った正弦波を作っておきます。

f:id:s0sem0y:20171024095433p:plain

scipyのsignal.lfilterに必要な係数の列 $a,b$ を以下のように渡してやれば、与えた差分方程式の計算を実行してくれます。これは、 $a$ が単なる $a[0]=1$ のみであるため、FIRフィルターなってます（最初にnp.convolveで見た奴と同じ）。

f:id:s0sem0y:20171024095527p:plain

実際、np.convolveで畳み込み計算を行うのと結果は全く同じになります（ただし、mode='full'では信号の両端を0パディングしている。右端は要らないデータであるため切り取っている）。

いろいろ試してみようと考えるのは良いのですが、以下の式

$\displaystyle \sum_{j=0}^{J-1} a[j]y[n-j] =\sum_{k=0}^{K-1} b[k]x[n-k]$

を睨んでみたところで、果たしてどうすればどんな処理が達成されるのか、というのを理解するのは相応に困難です。とりあえず、せっかくなので $a$ の方もしっかり系数列を与えてやることにしましょう。

f:id:s0sem0y:20171024100636p:plain

あら、ノイズが悪化しちゃいました。

んー、ほんの少しだけ値を変えてみましょう。

f:id:s0sem0y:20171024100843p:plain

さて、もともとの波形が跡形もなく飛んでいってしまいました。さて、なぜこのような事が起こったのかは分かるでしょうか。これは

$\displaystyle y[n] = \frac{1}{a[0]} \left( -\sum_{j=1}^{J-1} a[j]y[n-j] +\sum_{k=0}^{K-1} b[k]x[n-k] \right)$

この式の形式で見ると意外とすぐに分かります。 $y[n]$ はひとつ前の $y[n-1]$ に依存しているケースです。大雑把に書けば

$y[n] = -1.1 y[n-1] + ...$

という形式をしています。いいですか残りの入力がどんな値を持っていようが、前回の自分自身の出力の大きさを $1.1$ 倍するようなシステムですから、時間が立つに連れてとんでもない値になっていってしまうことが想像付くでしょう。

実はIIRフィルタは、設計を少し誤ると信号が発散する現象が起こってしまいます。一般にIIRフィルタには安定性が保証されていないという言い方をします。

FIRフィルタの方は大丈夫です。ココラヘンの議論をシステマチックに行うためには、更に話を進めていく必要がありますが、今回はこの辺で終わりにします。

最後に

周波数解析へ進むとできること

線形時不変システムを線形差分方程式で記述しておくことで、ある程度実現可能なシステムを統一的に扱うことが出来るようになりました。しかし、この段階ではまだまだフィルタの設計指針を掴んだとは言いづらいでしょう。

ここから更に進むためには、フーリエ変換・ラプラス変換・Z変換を通じた周波数解析が必要になってきます。周波数解析では、どのような係数の線形差分方程式を使うと、どのような周波数の振幅が強まるのか（あるいは弱まるのか）や位相がどれくらい遅れるか（あるいは進むか）などが定量的に評価されるようになります。

すなわち、信号をフーリエ変換した場合には、三角波の和として見ることになるわけですが、あるシステムが周波数 $ω$ の三角波にどのような影響を及ぼすのか、という形で議論するようになるわけです。システムが発散を起こすというのは、ある周波数の三角波の振幅が無限大になってしまうことで引き起こされます。

そして応用上重要な利点は、システムに、とある周波数に基づいて振幅の操作を行って欲しいと思った場合（要するにローパスフィルタだとかバンドパスフィルタを作りたい場合）、その周波数特性を先に決めてから、それを達成する線形差分方程式がどんな形になっているか（すなわち係数 $a,b$ がいくつか）を求める手段も手に入ります。

インパルス応答の議論は必要だったのか？

ところで、線形差分方程式で線形時不変システムの実用的なものをカバーできるならば、「インパルス応答との畳み込みが」という話は要らないのではないかと考えるかもしれません。

これは今、「線形時不変システムを設計しよう」という立場の説明が多いためだと思われます。仮に、未知の線形時不変システムがあり、それがどのような差分方程式になっているのか（あるいは微分方程式になっているのか）わからない時、そのシステムの振る舞いを知る手がかりとなるのがインパルス応答なのです。

仮にその未知のシステムが線形時不変システムであるならば（その微分方程式や差分方程式がどうであれ）、インパルス応答を見れば全ての情報が入っており、インパルス応答をしっかり把握できさえすれば、今後のシステムに対してどんな入力が来ても、その出力を予測する（つまり把握したインパルス応答との畳み込みを計算する）ことが可能です。

より具体的には、観測されたインパルス応答を達成するような線形差分方程式（線形微分方程式）を正確に設計するということが必要なケースがあるということです。もちろんインパルス応答を正しく計測できているのかなどの問題はあれど、手がかりになるのは間違いありません。

特に制御の問題では、システムが望んだ出力を達成するための入力を設計する必要があるため、未知のシステムがどのようなダイナミクスを持っているのかを把握すること（システム同定）は非常に重要です。システムの振る舞いが分かって初めて、では望みの出力をするような入力をどのように作るのかという話に移れるのです。

周波数領域での解析を最初からやればいいんじゃないのか

周波数領域での解析は非常に強力な設計指針を与えてくれます。実際、線形差分方程式を眺めていても分からないような情報が沢山得られます。

しかし、周波数領域での解析と時間領域での解析はどちらも重要です。

まず時間と周波数は、私達が解析してどのように理解するかは別として、物理的に結びついている量です。私達が周波数解析をしようがしまいが、周波数領域で何らかの事象を起こすためにはそれに対応する時間領域で事象があるはずです。それが時間領域でピタリと求まるならそれでも構いません。

特に音声処理の分野は、周波数スペクトル解析やケプストラム解析（周波数領域の周波数解析）などを行い、物理的に重要だと思われる特徴量の研究を進めてきましたが、リカレントネットワークは入力データの時間領域での振る舞いだけで上手く情報を持ってきてしまいますし（スペクトログラムを食わせるならば、ケプストラム解析に近いと思うけど）、wavenet(1D畳み込みをdilationしながら利用)なんかも凄まじい能力を持っています（ただ、マシンパワー無いと無理だから、ああ言うのは一般企業にはどっちしても発見できないと思いますが）。

あと状態空間表現というのは、時間領域で使われる非常に有用なモデルの一つです。

s0sem0y.hatenablog.com

HELLO CYBERNETICS

深層学習、機械学習、強化学習、信号処理、制御工学、量子計算などをテーマに扱っていきます

【信号処理・制御工学の基礎】線形差分方程式の表現まで

はじめに

前回の復習

インパルス応答でのシステムの考え方

インパルス応答を作る＝線形時不変システムを作る

補足

リアルタイムで動くシステム

線形差分方程式によるフィルタ設計

FIRフィルタとIIRフィルタ

FIRフィルタ

IIRフィルタ

線形差分方程式によるフィルタ処理

scipy.signal.lfilter

最後に

周波数解析へ進むとできること

インパルス応答の議論は必要だったのか？

周波数領域での解析を最初からやればいいんじゃないのか