dfilterh.dvi

12.5

デジタルフィルターの基礎

 離散時間信号を処理するシステム  線形時不変システム 信号値の乗算、加減算、及び信号値を記憶して遅らせるという3種類の処理の 組み合わせとして実現。ソフトウェア、ハードウェアのどちらでも実現でき、 多くの応用。  周波数特性の解析 システムの入力信号と出力信号との関係が、周波数に対してどのように変化す るかを調べること。  この資料は、参考文献：デジタル信号処理、貴家、昭晃堂、ISBN4-7856-1194-4 C3055からの引用。 12.5.1 信号の表現方法  正弦波信号 x(t)=Asin(t) は角周波数で、単位はrad/sec。  サンプリング間隔T s =1=F sでサンプリング。 t=nT sにおける「 n番目」の 信号値は、 x(nT s )=Asin(nT s ) ただし、F sはサンプリング周波数。  正規化表現 !=T sとおき、 x(nT s )!x(n)と表すと、 x(n)=Asin(nT s )=Asin(! n) ここで!は正規化角周波数で、単位はrad。  正規化周波数 { 周波数F(Hz)に対して、=2F { 正規化周波数f を、!=2fとして定義。f の単位は無次元。 { !と、fとF の間の関係 ! = T s ==F s (1) f = !=(2 )=F=F s (2) 12.5.2 移動平均フィルターの例  図2は図1の連続信号をサンプリングしたもの。これを左からx(0);x(1);x(2);::: と表す。  移動平均システム 図3は、図2のx(0);x(1);x(2);:::に対して、一つ前の信号値との平均を取っ て「作った」信号。つまり、 y (n)= 1 2 (x(n)+x(n01)) 図3の波形は、図2よりも高周波成分が小さい。

 図4は、連続する3点の移動平均 y (n)= 1 3 (x(n)+x(n01)+x(n02)) 図3よりもさらに高周波成分が小さい。  ローパスフィルター 移動平均をとると、高周波成分が小さくなり、低周波成分だけが残る。つまり、 信号x(n)→[移動平均]→低周波成分y(n)  デジタルフィルター 一般に、サンプリング(と量子化)がなされた離散信号を入力として、演算処理 を行い、新しい離散信号を出力するものをデジタルフィルターという。 図 1: 元の信号 図 2: サンプリングされた信号 図3: 2点平均化された信号 図 4: 3点平均化された信号  正弦波に対する移動平均 図5、6は、正弦波信号と、その5点平均と7点平均とを重ね書きしたもの。こ れらから分かることは、 { 出力波形も同じ周期(周波数)正弦波である。 { ただし、振幅が少し小さく、また、波が全体に右に寄っている(位相が遅 れている)。  以上の二つの例から、 { なぜ平均処理により信号の大きさが変わり、位相がずれるのか? { 平均をとる点数と、大きさや位相の変動との関係は?

図 5: 正弦波と5点平均 図 6: 正弦波と7点平均 12.5.3 準備：基本的な信号とその性質  正弦波信号 sin(! n)または cos(!n)  複素正弦波信号 e i! n =cos(! n)+isin(! n) i は虚数単位 1 。この関係はオイラーの公式と呼ばれる。 この複素正弦波信号は、後述する「周波数特性」を扱う上で重要。 これを用いて、（実数の）正弦波信号を複素正弦波信号で表すことができる。 cos(! n) = e i! n +e 0i! n 2 sin(!n) = e i! n 0e 0i! n 2 (3)  単位ステップ信号u(n)(図7(a)) u(n)= ( 1 n0 0 n<0  インパルス(信号)(n)(図7(b)) (n)= ( 1 n=0 0 n6=0 図7: 単位ステップとインパルス

図 8: インパルス(n)の性質 12.5.4 インパルスの性質  例)(n02)という信号はn=2の時だけ1になる信号(図8の(a))  図8の(b)のように、x(01)=01; x(0)=2; x(2)=1でそれ以外は0の信 号x(n)は、インパルスを用いて次のように書くことができる。 x(n)= 上の式は、 x(n)=x(01) (n+1)+x(0)(n)+x(1)(n01)+x(2) (n02) と書くこともできる。  一般に、信号x(n)は、 x(n)= 1 X k=01 x(k) (n0k ) (4) のようにインパルスを用いて級数の形で表現できる。 12.5.5 線形・時不変システム  次に信号の変換について考えよう。信号処理システムを入力信号x(n)を出力 信号y (n)に一意的に変換するものとして定義し、 y(n)=T[x(n)] と表す。 変換の条件によって、システムを以下のように分類できる。  時不変(シフト不変)システム(time-invariantsystem) 入力x(n)に対応する出力をy (n)とするとき、任意の整数kについて y(n0k)=T[x(n0k )] が成立するようなシステム。  例）図9(a)(b) 1 工学系では記号としてjを使うことが多い。

図9: 線形・時不変システムの性質  線形システム(linear system) 任意の入力x 1 (n);x 2 (n)に対応する出力をそれぞれy 1 (n);y 2 (n)とするとき、 任意の定数a;bに対して T[ax 1 (n)+bx 2 (n)] = aT[x 1 (n)]+bT[x 2 (n)] = ay 1 (n)+by 2 (n) が成立するようなもの。  例)図9(c)  線形時不変システム 線形性と時不変性の条件を同時に満たすとき、線形時不変システム(linear time-invariantsystem)と呼ばれる。線形性と時不変は、独立な条件。  因果性システム(causalsystem) 任意の時刻n 0における出力 y(n 0 )が、その時刻よりも過去の時刻n≦n 0のみ の入力x(n)を用いて計算されるシステムである。 例)因果性システムではないもの

12.5.6 線形時不変システムにおけるたたみ込みとインパルス応答  インパルス応答(impulseresp onse)

システムにインパルスを入力した場合の出力h(n)： h(n)=T[ (n)]  線形時不変システムでは、システムのインパルス応答h(n)がわかると、任意 の入力x(n)に対する出力を次式で求めることができる(証明は後述)。 y(n)= 1 X k=01 x(k )h(n0k) (5) あるいは、m=n0kとおくと、 y (n)= 01 X m=1 x(n0m)h(m)= 1 X k =01 h(k )x(n0k ) (6)  たたみ込みによるシステムの表現 上式の関係を信号x(n)とh(n)のたたみ込み(convolution) といい、次のよう に書く。 y(n)=x(n)?h(n)  例）3点平均とたたみ込み インパルス応答がh(0)=h(1)=h(2)= 1 3 のシステムを求めて見よう。  たたみ込みの導出 式(5)のたたみ込みの導出のためには、線形性と時不変性の条件が必要である。 まず、システムの入出力関係は、y (n)=T[x(n)]と表現される。これに式(4) を代入する。 y(n)=T[x(n)]=T[ 1 X k =01 x(k )(n0k)] 線形性を仮定すると、上式は y (n) = 1 X k =01 T[x(k )(n0k )] (7) = 1 X k =01 x(k )T[(n0k )] (8) 式(7)から式(8)の変形で、 P の中でx(k)は定数であることに注意。 インパルス応答をh(n)=T[ (n)]とすると、時不変性から、 h(n0k )=T[ (n0k )] これを式(8)に代入して、式(5)のたたみ込みの式を得る。  例) 図10

図 10: たたみ込みの例 12.5.7 ハードウェアでの実現  演算要素 式(5)からわかるように、たたみ込みは乗算、加減算、シフト(遅延)という三 つの要素によって実現できる。  図11の上図はこれらの要素を示し、図11の下図は、次のシステムを実現した もの。 y (n)=x(n)+2x(n02) 12.5.8 再帰型システム  次のシステムを考えよう。 y(n)=x(n)+by(n01) これは右辺にyが現れており、たたみ込みの形をしていない。しかし、このシ ステムをハードウェアで実現することは可能で、それは図12のようになる。こ こでは、ある時刻の出力を入力として用いてその後の出力を計算している。こ のような処理をフィードバックといい、フィードバックを含むシステムを再帰 型システムという。  インパルス応答 この図において信号の流れを追いかけることによって、インパルス応答を求め てみると、時刻n=0から始めて、出力は、1;b;b 2 ;b 3 ;111と出力が無限に続く ことが分かる。  I IR とFIR

このように無限個のインパルス応答をもつシステムをIIR (innite impulse resp onse) システムといい、逆に有限個のインパルス応答を持つシステムを FIR(niteimpulseresp onse)システムという。移動平均などのシステムはFIR

であった。I IRシステムは再帰型であるが、再帰型であってもI IRとならない

ものもある。

 前述のようにインパルス応答が分かると、たたみ込みの式(6)によってシステム (入力x(n)に対する出力y (n))を記述することができた。再帰型のシステムで

図 11: ハードウェアによる実現

あってもそのようにできる。このシステムのインパルス応答はh(k)=b k ;(k= 0;1;2;:::)であるから、システムは、 y(n)= 1 X k =0 b k x(n0k) と無限級数で表現できる。  I IRシステムでは安定性に注意 安定性とは、有限の信号値をもつ信号をシステムに入力したときに、出力が必 ず有限の値となるようなシステムである。図12のシステムでは、もしjbj>1 であると、出力が無限大になり、不安定なシステムとなる。 12.5.9 システムの周波数特性 システムを利用するためには、その特性について調べておく必要がある。このよう な解析を行うには、信号の周波数に着目し、周波数に対するシステムの特性を求めて やる。これを調べるには、システムに対して、正規化角周波数!の(複素)正弦波信 号を入力し、その出力を求めればよい。 そして、もしすべての周波数についてシステムの周波数特性を調べることができれ ば、任意の正弦波信号に対する出力を容易に知ることが可能となる。また、詳細は次 節で述べるが、正弦波信号以外の入力に対する出力も、この周波数特性から決定され る。つまり、インパルス応答と同様に、周波数特性は、線形時不変システムを完全に 記述している。  振幅特性と位相特性 線形時不変システムでは、正弦波信号x(n)=cos(!n)を入力したとき、出力 信号は、

y(n)=A(! )cos(!n+ (!))

すなわち次のような特徴がある。 { 出力も同じ周波数(!)を持つ正弦波信号である。 { システムは正弦波信号の大きさ(A(! ))と位相((! ))のみを変える働きが ある。 { 大きさ(A(!))と位相( (! ))は、周波数(!)の関数であり、入力信号の周 波数によって決まる。 システムの周波数特性(frequencycharacteristics)とは、次の振幅特性と位相 特性の両方で表現される。

{ 入力と出力の大きさ(振幅)の関係A(!)を振幅特性(amplitude character-istics)といい、 { 位相の関係(! )を位相特性(phasecharacteristics)という。  複素正弦波信号入力 次に、複素正弦波信号 e i! n =cos(! n)+isin(! n) に対する出力を考えてみよう。システムのインパルス応答h(n)は与えられて いると仮定すると、出力を求めるにはたたみ込みの式(6)に上式を代入すれば

よい。整理すると、 y (n) = 1 X k =01 h(k )e i! (n0k ) = e i! n 1 X k=01 h(k)e 0i! k (9) = e i! n H(e i! ) (10) ここで、Hは、式(9)の6の部分は、 H(e i! )= 1 X k =01 h(k)e 0i! k (11) H をe i! の関数と捉えていることに注意。 このHの値は複素数であるので、 H(e i! )=A(! )e i (! ) (12) のように、大きさA(!)と偏角 (!)を用いて極座標表現できる。この式を式 (10)に代入すると、 y (n) = e i! n H(e i! )=A(! )e i(! n+ (! )) つまり、 { 複素正弦波信号に対する出力信号は、振幅と位相のみが入力信号と異なる。 { インパルス応答が既知であれば、式(11)を計算し、式(12)へ変形するこ とによって、大きさと位相の特性を知ることができる。（位相の特性を知 るために入力信号は不要）  例）3点移動平均y (n)= 1 3 (x(n)+x(n01)+x(n02))の振幅特性と位相特 性は、 図13: 3点移動平均の周波数特性

位相特性のグラフがノコギリ歯のようになっているのは、e i! =e i(! +2 ) が成立 することから、0 (! )<の範囲で位相特性(! )を描いているからであ る。ノコギリ歯のようにせず、01 (!)<1の範囲で直線で描いてもよい。  正弦波信号に対する周波数応答との関係 複素正弦波信号は複素数であり、現実的に我々は正弦波信号を与えたときの周 波数特性に興味がある。システムに正弦波信号を与えたときの周波数特性は、 複素正弦波信号の周波数特性と基本的に一致する。  3点移動平均の周波数応答 図13は、正弦波信号に対する応答も表していることから、次のことが言える。 { !=の付近でA(!)の絶対値が小さくなっている。!=は、周波数で 言うと、F= ! 2 p i F sに !=を代入してF =F s =2となる。つまり、サン プリング周波数の半分あたりで、振幅が小さくなることを示している。こ れは、先に述べたように高い周波数成分をカットするローパスフィルター の働きを示している。 { 位相特性は、周波数に対して線形であり、周波数が高いほど位相が遅れる ことを示している。 12.5.10 z変換と伝達関数 ここまでで見たように、システムの周波数特性を求めるには、H(e i! )を求めれば よい。前項では、その一つの方法としては、システムのインパルス応答を求めて、式 (11)を計算し、それから式(12)を導出した。しかし、フィードバックを含む再帰型 システムでは、I IRシステムとなることがあり、インパルス応答は無限の長さにな り、面倒である。そこで通常は、以下で紹介するようなz変換という方法を用いて、 H(e i! )を求めている。 本節では、z変換と、この変換法に基づくシステムの伝達関数とシステムの周波数 特性を求める方法について説明する。  離散時間信号x(n)のz変換X(z)を次式で定義する。 X(z)= 1 X n=01 x(n)z 0n (13) ここで、zは複素変数である。いま、表現を簡潔にするために、x(n)のz変換 がX(z)であるとき両者の関係を x(n) z $X(z) または Z[x(n)]=X(z) と表す。下では、時間信号を小文字(x(n))、変換された信号を大文字(X(z ))で 表現する。  例）インパルス(n)のz変換 図14(a)のインパルス(n)のz変換は1。

 例）図14(b)の信号の2(n02)のz変換 これから、一般にc(n0k )のz変換はcz 0k となることが分かる。ただし、c は定数。  例)図14(c)の信号x(n)のz変換 任意の信号x(n)はインパルスを用いて表現できることから(式(4))、この例の ように、インパルスに対するz変換から容易に任意の信号のz変換を求めるこ とができる。 図14: 信号例  z変換の性質 z変換を実際の場面で使いこなすためには、以下に示すz変換の性質を理解す る必要がある。z変換では、常に以下の関係が成立する。 { 線形性 任意の2つの信号x 1 (n);x 2 (n)のz変換をそれぞれX 1 (z)=Z[x 1 (n)];X 2 (z)= Z[x 2 (n)]とする。このとき、 Z[ax 1 (n)+bx 2 (n)] = aZ[x 1 (n)]+bZ[x 2 (n)] = aX 1 (z)+bX 2 (z) が成立する。ここで、aおよびbは任意の定数である。 { 時間シフト 信号x(n)のz変換がX(z )=Z[x(n)]であるとき、 Z[x(n0k)]=X(z)z 0k ; (k:任意の整数)

が成立する。  たたみ込みとz変換の積 任意の2つの信号x 1 (n);x 2 (n)のz変換をそれぞれX 1 (z);X 2 (z)とする。この とき、両者がたたみ込みの関係にあるとき(式(5))、 x 1 (n)?x 2 (n)= 1 X n=01 x 1 (n)x 2 (n0k) z $X 1 (z)X 2 (z) (14) が成立する。これは次のようにして導くことが出来る。 y (n)= 1 X k =01 x 1 (k )x 2 (n0k )とおいて、これをz変換すると、 Z[y (n)] = 1 X n=01 y(n)z 0n = 1 X n=01 ( 1 X k=01 x 1 (k )x 2 (n0k ))z 0n = 1 X k =01 x 1 (k )( 1 X n=01 x 2 (n0k )z 0n ) ここで、n0k =pとおくと、z 0n =z 0p z 0k となり、これを代入して整理す ると、 = 1 X k =01 x 1 (k)( 1 X p=01 x 2 (p)z 0p z 0k ) = 1 X k =01 x 1 (k)z 0k ( 1 X p=01 x 2 (p)z 0p ) = 1 X k =01 x 1 (k)z 0k X 2 (z) = X 1 (z)X 2 (z )  システムの伝達関数 ここでは、インパルス応答に代わる表現として、システムの伝達関数を定義す る。I IRシステムのように、インパルス応答では無限の表現が必要になるシス テムに対しても、伝達関数は有限な表現を与え、表現としてより簡潔である。  伝達関数の定義 システムの伝達関数をシステムのインパルス応答h(n)のz変換H(z)として定 義しよう。すなわち、 H(z )=Z[h(n)] さて、線形時不変システムでは、入力信号x(n)、インパルス応答h(n)と出力 信号y (n)の間には、たたみ込みの関係が成立した。 y(n)=x(n)?h(n)= 1 X k =01 h(k )x(n0k ) z変換のたたみ込みの性質から、これをz変換すると、 Y(z)=H(z )X(z ) と表現できる。ただし、Y(z) =Z[y (n)];H(z )=Z[h(n)];X(z)=Z[x(n)]で ある。ここで、H(z)をシステムの伝達関数という。

したがって、伝達関数H(z)は、入出力信号のz変換の比： H(z)= Y(z) X(z) としても計算できる。  非再帰型システムの伝達関数 まず、例として、3点平均を求めるシステムで伝達関数を二つの方法で求めて 見よう。 y (n)= 1 3 (x(n)+x(n01)+x(n02)) のインパルス応答はh(n)= 1 3 ((n)+(n01)+(n02))であるので、伝達関 数は、h(n)のz変換として、 H(z)= 1 3 (1+z 01 +z 02 ) 第2の方法として入出力信号のz変換の比として求めるには、まずZ[y (n)]= Y(z)、Z[x(n)] = X(z)とする。さらに、時間シフト性よりZ[x(n01)] = X(z)z 01 ;Z[x(n02)]=X(z)z 02 だから、 Y(z) = 1 3 (X(z)+X(z)z 01 +X(z)z 02 ) = 1 3 (1+z 01 +z 02 )X(z ) これから、 H(z)= Y(z) X(z) = 1 3 (1+z 01 +z 02 ) 両者のアプローチに本質的な差異はない。しかし後述のように再帰型システム に対しては、後者がより簡潔に伝達関数を与える。  伝達関数の一般形 たたみ込みの式(6)から、因果性を満たす非再帰型システムは一般に、次のよ うに与えられる2 y (n)= N01 X k =0 h(k )x(n0k ); (N :正整数) 3点平均の場合と同様に、伝達関数を求めると、 Y(z)= N01 X k=0 h(k )X(z)z 0k =( N01 X k=0 h(k )z 0k )X(z) よって、 H(z)= Y(z) X(z) = N01 X k =0 h(k )z 0k 伝達関数のzの係数が直接インパルス応答に対応することに注意してほしい。ま た、伝達関数はzの多項式であり、この多項式の次数を伝達関数の次数(order) という。この伝達関数の次数はN01次である。  z変換による応答の計算 伝達関数H(z )のシステムに信号X(z)を入力したときの出力はH(z)X(z)に なることを利用して、x(n)に対する出力を計算することができる。 2 k<0でインパルスは信号値が0なので因果性より応答も0、つまりk<0ではh(k)=0。また非 再帰的であるからFIRとなり、ある時刻N>0より先の応答は0。つまりk>Nでh(k )=0。

 例)x(0)=2; x(1)=02; x(2)=2、それ以外でx(n)=0なる入力を2点平均 システムy (n)= 1 2 (x(n)+x(n01))に入力したときの出力を求めると、 このようにして、z変換を用いて出力を計算することができる。より一般には 逆z変換を用いる必要がある。  再帰型システムの伝達関数の導出法 フィードバックを持つ再帰型システムとして、 y(n)=x(n)+by(n01) を例にする。先の非再帰型システムの場合とほぼ同様に、システムの伝達関数 を求めることができる。まず、両辺をz変換すると Y(z)=X(z )+bY(z)z 01 となる。次に、Y(z)とX(z)について整理すると Y(z)(10bz 01 )=X(z) を得る。ゆえに、伝達関数H(z)は、 H(z)= Y(z) X(z) = 1 10bz 01  伝達関数と周波数特性 周波数特性を伝達関数から計算することができる。伝達関数H(z)が既知であ るとき、そのzにe i! を代入する。すなわち、 H(e i! )=H(z)j z =e i! この方法の正当性は、次のように説明される。まず、入力とする信号をx(n)と して、それをz変換(式(13))する。 X(z)= 1 X n=01 x(n)z 0n 定義よりインパルス応答h(n)のz変換が伝達関数H(z)であるので、このx(n) にh(n)を代入すればH(z)が求まる。 H(z)= 1 X n=01 h(n)z 0n ここで、このzにe i! を代入すると、 H(e i! )= 1 X n=01 h(n)e 0i! n となり、式(11)と一致する。すなわち、伝達関数のzにe i! を代入した結果と、 直接にインパルス応答から式(11)を計算した結果とは一致する。 Z変換による周波数応答についてまとめたものを図15に示す。

ᤨ㑆㗔ၞ

y(n)=T[x(n)]

๟ᵄᢙ㗔ၞ

๟ᵄᢙᔕ╵

H(exp(iȁ))

\㗔ၞ

વ㆐㑐ᢙ

H(z) = Y[z]/X[z]

ࠗࡦࡄ࡞ࠬᔕ╵

h(n)

y(n) = Ǜh(k)x(n-k)

x(n) = exp(iȁn)

zᄌ឵

z = exp(iȁ))

図15: z変換と周波数応答 12.5.11 まとめ  デジタルフィルター 簡単な演算によってローパスフィルターなどの信号処理を行うことができる。  周波数解析 デジタルフィルターの周波数特性を求める方法  デジタルフィルターの設計 所望の周波数特性をもつデジタルフィルターを合成すること。

図 16: 各種のフィルター（LPF (Low Pass Filter), BPF(Band Pass Filter), HPF (HighPassFilter),BRF(Band RejectFilter)

プログラム例

// lpf1k.java // 低域濾過フィルター // エンコーディングは次の仕様に固定 // サンプリングレート8kHz、データ長 8ビット、モノラル、PCM // (C) 鈴木宏正＠東京大学 1999/9/13 略

static void convert(String infn, String outfn) { FileInputStream fin; FileOutputStream fout; //低域通過タイプ //減衰量: 60.000 [dB] //次数: 100 次 //標本化周波数: 8.00000 [kHz] //遮断周波数１: 1.00000 [kHz] double h[] ={ 0.00012980, 0.00012571, 0.00000000,-0.00021187,-0.00037517, -0.00032645, 0.00000000, 0.00047473, 0.00079661, 0.00066250, 0.00000000,-0.00089617,-0.00146002,-0.00118280, 0.00000000, 0.00153030, 0.00244607, 0.00194765, 0.00000000,-0.00244520, -0.00385770,-0.00303525, 0.00000000, 0.00373275, 0.00583714, 0.00455654, 0.00000000,-0.00553114,-0.00860537,-0.00668972, 0.00000000, 0.00807882, 0.01255811, 0.00976658, 0.00000000, -0.01185753,-0.01853088,-0.01452100, 0.00000000, 0.01805053, 0.02870358, 0.02300258, 0.00000000,-0.03057276,-0.05112438, -0.04387509, 0.00000000, 0.07433807, 0.15850612, 0.22485221, 0.25000417, 0.22485221, 0.15850612, 0.07433807, 0.00000000, -0.04387509,-0.05112438,-0.03057276, 0.00000000, 0.02300258, 0.02870358, 0.01805053, 0.00000000,-0.01452100,-0.01853088, -0.01185753, 0.00000000, 0.00976658, 0.01255811, 0.00807882, 0.00000000,-0.00668972,-0.00860537,-0.00553114, 0.00000000, 0.00455654, 0.00583714, 0.00373275, 0.00000000,-0.00303525, -0.00385770,-0.00244520, 0.00000000, 0.00194765, 0.00244607, 0.00153030, 0.00000000,-0.00118280,-0.00146002,-0.00089617, 0.00000000, 0.00066250, 0.00079661, 0.00047473, 0.00000000, -0.00032645,-0.00037517,-0.00021187, 0.00000000, 0.00012571, 0.00012980 }; 略

for (i = 0; i < size; i++) { byte temp; buf[0] = (byte)fin.read(); o = 0.0; for (m = 0; m<=NF; m++) { o += h[m]*((double)buf[m]); } fout.write((byte)o); for (m=NF; m>0; m--) { buf[m] = buf[m-1]; } } 略