趣味人のブログ

Cogito, ergo sum. 我思う故に我あり.

13. FIR フィルタの概要

本章では,13.1 において FIR (Finite Impulse Responce) フィルタ (有限インパルス応答フィルタ,非巡回型フィルタ) の基本を,13.2 において基本的な FIR フィルタの例を,13.3 において窓関数による FIR フィルタの例を述べる.次に,13.4 において,位相歪が生じない線形位相フィルタの直感的な原理を述べ,13.5 において線形位相フィルタのインパルス応答の例を示す.そして,13.6 においてこれに基づいた線形位相 FIR フィルタの例を,13.7 において窓関数による線形位相 FIR フィルタの例を示す.これらは教科書に十分な説明が記載されているため,本資料は留意すべき点を Note に示す.

13.1. FIR フィルタの基本

システム伝達関数 H (z) が以下の式で与えられるディジタルフィルタの構造を示す.ここで,h [n], (n ≥ 0) はディジタルフィルタのインパルス応答を,X (z), Y (z) は各々ディジタルフィルタの入出力となる離散時間信号 x [n], y [n], (n ≥ 0) の Z 変換を示す.

(13-1) Hz=YzXz=n=0N-1hnz-n

上記のディジタルフィルタに入出力される離散時間信号 x [n], y [n] の関係は,以下の畳み込みで記述される.

(13-2) yn=k=0N-1hkxn-k

FIRフィルタとは,上記の畳み込みによる離散時間信号 x [n], y [n] の関係をフィルタの実現アルゴリズムと見做したものであり,出力の離散時間信号 y [n] を現在及び過去 N − 1 回の入力 x [n], ..., x [n − (N − 1)] から求める.

FIR フィルタの基本構造を図13‑1に示す.図において,四角は入力信号に 1 サンプリング周期の遅延を加えて出力するレジスタ,三角は定数係数の乗算器,丸は加算器を示し,上記の畳み込みの関係を式の通りに計算する.レジスタの初期値は全て 0 となる.ここで,FIR フィルタの次数 N − 1 は,レジスタ数となる.

f:id:suzumushi0:20170611122207p:plain
図13-1: FIR フィルタの基本構造.

上記の基本構造から明らかな通り,FIR フィルタの設計にはインパルス応答を求める必要がある.フィルタが満たすべき周波数伝達関数Hd (e jωT) とすると,式 (12‑11) に示した離散時間フーリエ変換 (DTFT) の逆変換からインパルス応答を求める事が出来る.但し,式 (12‑11) の積分を求める事が容易では無い場合は,周波数伝達関数をサンプリングして,式 (12‑4) に示した離散フーリエ変換 (DFT) の逆変換からインパルス応答の近似を求めても良い.勿論,フィルタが満たすべきシステム伝達関数の逆 Z 変換からインパルス応答を求めても良い.

多くの場合,フィルタが満たすべき周波数伝達関数 Hd (e jωT) に対応したインパルス応答 h [n] は無限に続くため,Hd (e jωT) と h [n] の関係は以下で記述される.

(13-3) HdejωT=n=0hne-jnωT

ディジタル信号処理では,入出の離散時間信号 x [n] が出力の離散時間信号 y [m], (n > m) に影響を及ぼす事はあり得ないため,インパルス応答は必ず h [n] = 0, (n < 0) となる.また,ディジタル信号処理では,初期条件を 0 とする事を暗黙の前提としている.従って,上記の周波数伝達関数は式 (12‑10) に示した離散時間フーリエ変換 (DTFT) の定義式と等価となる

一方,FIR フィルタはインパルス応答 h [n] を h [N − 1] で打ち切るため,FIR フィルタによって得られる現実の周波数伝達関数 H (e jωT) は以下で与えらえる.

(13-4) HejωT=n=0N-1hne-jnωT

ここで,インパルス応答を打ち切る事は,離散時間フーリエ変換 (DTFT) の総和 Hd (e jωT) を部分和 H (e jωT) で近似する事と等価となる.12.2 に述べた通り,フーリエ級数と離散時間フーリエ変換 (DTFT) は時間と角周波数と入れ替えた関係にあるから,この近似によって FIR フィルタの周波数特性には,8.1 に述べたギブスの現象に起因するリップルが生じる

13.2. 基本的な FIR フィルタの例

11.2 の例と同じインパルス応答 h [n] に基づいた FIR フィルタの例を示す.これは,時定数 τ のローパスフィルタのインパルス応答を周期 T でサンプリングしたものである.式 (11‑8) に示したインパルス応答 h [n] を以下に再掲する.

(13-5) hn=1τe-nTτ,n0

上記のインパルス応答に基づいた FIR フィルタのシステム伝達関数 H (z) を以下に示す.カットオフ周波数は ωc = 1 / τ とする.また,FIR フィルタの次数 (レジスタ数) はN − 1 とする.

(13-6) Hz=n=0N-1hnz-n=1τn=0N-1e-nTτz-n=ωcn=0N-1e-nωcTz-n

上記のシステム伝達関数の直流におけるゲインを 1 倍とした正規化システム伝達関数HN (z) 以下で与えられる.(式の展開は省略,等比級数の和の公式を使用,以下同様).

(13-7) HNz=Hzlimz1Hz=1-e-ωcT1-e-NωcTn=0N-1e-nωcTz-n

上記の正規化システム伝達関数に対応する,周波数伝達関数を以下に示す.(式の展開は省略).

(13-8) HNejωT =1-e-ωcT1-e-NωcTn=0N-1e-nωcTe-jnωT =1-e-ωcT1-e-NωcT1-e-NωcT+jωT1-e-ωcT+jωT

N = 16 の場合の上記の周波数伝達関数のボード線図を図13‑2に示す.横軸の角周波数は ω / ωc として正規化されている.図中の青線は正規化されたナイキスト周波数 fn を示しており,この例ではカットオフ周波数 ωc の 20 倍の角周波数としている.また,図中の赤線はアナログフィルタの特性を示す.

f:id:suzumushi0:20170611122208p:plain
図13-2: FIR によるローパスフィルタのボード線図の例.

図から明らかな通り,ゲイン線図,位相線図共に,インパルス応答の打ち切りによって生じるギブスの現象のため,周波数特性が波打つリップルが生じている.

13.3. 窓関数による FIR フィルタの例

8.1 に示した通り,窓関数によってギブスの現象に起因するリップルを抑制する事が出来るため,上記の FIR フィルタにハミング窓を適用した例を示す.この例では n = 0 周辺の項の重み付を高くする必要があるため,式 (8‑17) で定義した以下のハミング窓を適用する.

(13-9) wn=0.54+0.46cosπnN-1,0nN-1

この場合のインパルス応答 h [n] は以下で与えられる.

(13-10) hn=1τ0.54+0.46cosπnN-1e-nTτ,n0

上記のハミング窓を適用したインパルス応答に基づいた FIR フィルタのシステム伝達関数 H (z) を以下に示す.(式の展開は省略).

(13-11) Hz =n=0N-1hnz-n =ωc0.54n=0N-1e-nωcTz-n +0.23n=0N-1e-nωcT-jπN-1+e-nωcT+jπN-1z-n

上記のシステム伝達関数の直流におけるゲインを 1 倍とした正規化システム伝達関数は以下で与えられる.(式の展開は省略).

(13-12) HNz=Hzlimz1Hz limz1Hz=ωc0.541-e-NωcT1-e-ωcT +0.231-e-NωcT-jπN-11-e-ωcT-jπN-1+1-e-NωcT+jπN-11-e-ωcT+jπN-1

上記の正規化システム伝達関数に対応する周波数伝達関数を以下に示す.(式の展開は省略).

(13-13) HNejωT =ωclimz1Hz0.541-e-NωcT+jωT1-e-ωcT+jωT +0.231-e-NωcT-jπN-1+jωT1-e-ωcT-jπN-1+jωT +1-e-NωcT+jπN-1+jωT1-e-ωcT+jπN-1+jωT

上記の周波数伝達関数のボード線図を図13‑3に示す.図13‑2と同様に N = 16 であり,横軸の角周波数は ω / ωc として正規化されている.また図中の青線は正規化されたナイキスト周波数 fn を示しており,カットオフ周波数 ωc の 20 倍の角周波数としている.

f:id:suzumushi0:20170611122209p:plain
図13-3: 窓関数による FIR フィルタのボード線図の例.

図から明らかな通り,ゲイン線図,位相線図共に,ハミング窓によって周波数特性のリップルが大幅に抑制されている.その反面,ハミング窓による重み付の影響で,減衰量が低下している.

以上より,IIR フィルタで出来る事を FIR フィルタで実現する事はあまりお勧めできない.11.2 に述べた IIR によるローパスフィルタに必要なレジスタ数は,式 (11‑11) の差分方程式から明らかな通り,僅か 1 個である.一方,上記の FIR フィルタの次数,即ちレジスタ数は N − 1 = 15 個となる.これはフィルタをハードウエアで実現する場合は論理回路のリソース数,ソフトウエアで実現する場合は実行時間に大幅な差が生じる事を示している.また,上記から明らかな通り FIR フィルタの設計は IIR フィルタと比較して煩雑である.従って,FIR フィルタが本領を発揮するのは,次に述べる線形位相 FIR フィルタの様な IIR フィルタでは実現できない特性のフィルタとなる.

13.4. 線形位相 FIR フィルタの基本

8.2 に述べた通り,正弦波を合成した際に得られる波形は,各正弦波の振幅と初期位相に依存する.従って,フィルタのゲイン特性が角周波数によらず 1 倍であっても,信号の遅延時間が角周波数によって異なる場合は,入力信号と異なる波形の出力信号が得られる.これを位相歪と呼ぶ.理想的にはフィルタの位相特性は角周波数によらず 0,即ち位相の進み遅れが無い零位相特性となる事が望ましい.しかし,この様なフィルタを実現する事は出来ないため,位相歪を避ける必要がある場合は,フィルタの遅延時間が角周波数によらず一定となる線形位相 (直線位相) 特性のフィルタが用いられる.

零位相特性フィルタ

位相特性が角周波数によらず 0 となる零位相特性フィルタの直感的な原理を説明する.インパルス応答を h [n] (n ≥ 0),離散時間信号 x [n], y [n], (n ≥ 0) を入出力とするディジタルフィルタを考える.この場合,離散時間信号 x [n], y [n] の関係は,以下の畳み込みで記述される.また,入出力の離散時間信号 x [n], y [n] の周波数成分の間には,周波数伝達関数の位相特性に応じた位相の遅れ進みが生じる.

(13-14) yn=k=0hkxn-k

ここで,入力の離散時間信号 x [n] の時間を逆順にした離散時間信号 rx [n] = x [−n], (0 ≥ n) を考える.この離散時間信号 rx [n] を上記のフィルタの入力とし,出力の離散時間信号を ry [n], (0 ≥ n) とすると,rx [n], ry [n] の関係は,以下の畳み込みで記述される.また,入出力の離散時間信号 rx [n], ry [n] の周波数成分の間には,上記と同じ位相の遅れ進みが生じる.

(13-15) ryn=k=0hkrxn-k=k=0hkx-n+k

次に,これらの離散時間信号 rx [n], ry [n] の時間を逆順にした離散時間信号 rrx [n] = rx [−n] 及び rry [n] = ry [−n], (n ≥ 0) を考える.ここで rrx [n] = x [n] は明らかである.よって,離散時間信号 x [n] と rry [n] の周波数成分の間には,上記と逆の位相の遅れ進みが生じる従って,入力の離散時間信号 x [n] と y [n] + rry [n] の周波数成分の間には位相の遅れや進みは生じない

ここで,h [n] の時間を逆順とした理論上,想像上のインパルス応答を rh [n] = h [−n], (0 ≥ n) とすると,上記の離散時間信号 x [n], rry [n] の関係は,以下の畳み込みで記述される.

(13-16) rryn=ry-n=k=0hkxn+k=k=-0rhkxn-k

上記の式から明らかな通り,散時間信号 rry [n] は,入力の離散時間信号 x [n] と rh [n] = h [−n] の畳み込みから得られる.よって,入力の離散時間信号を x [n],インパルス応答を h [n] + rh [n] = h [n] + h [−n] とするフィルタの出力の離散時間信号は y [n] + rry [n] となる

従って,入力の離散時間信号を x [n], (n ≥ 0),出力の離散時間信号を y [n],理論上,想像上のインパルス応答を hd [n] = hd [−n] とするフィルタの入出力の関係は以下の畳み込みで記述され,この入出力の離散時間信号間には位相の進みや遅れは生じず,即ち位相特性は角周波数によらず 0 となる.このため,想像上,理論上のインパルス応答が hd [n] = hd [−n] となるフィルタの位相特性は零位相特性となる

(13-17) yn=k=-hdkxn-k

線形位相 (直線位相) FIR フィルタ

残念ながら,13.1 に述べた通り,ディジタル信号処理のインパルス応答は h [n] = 0, (n < 0) でなけれならないため,零位相特性フィルタを実現する事は出来ない.しかし,FIR フィルタでは,信号の遅延時間が角周波数によらず一定となり,位相歪が生じない線形位相 (直線位相) 特性を実現する事ができる.零位相特性となるインパルス応答 hd [n] に対して,以下のインパルス応答 h [n] を考える.ここで,N は奇数とする.

(13-18) hn= hdn-N-12,0nN-10,n<0,nN

上記から明らかな通り,h [n] は,hd [n] を |n| ≤ (N − 1) / 2 で打ち切り,かつ (N − 1) / 2 サンプル時間遅延させたインパルス応答となる.ここで,hd [n] = hd [−n] であるから,h [n] には以下の関係が成立する.

(13-19) hn=hN-1-n,0nN-1

従って,インパルス応答 h [n] が式 (13‑19) を満たす FIR フィルタの遅延時間は角周波数よらず一定となり,線形位相特性となる.線形位相 FIR フィルタの構造を図13‑4 に示す.

f:id:suzumushi0:20170611122210p:plain
図13-4: 線形位相 FIR フィルタの構造.

Note: N が偶数の場合でも,インパルス応答が式 (13‑19) を満たす FIR フィルタは線形位相特性となる.尚,この場合の FIR フィルタの構造は図13‑4 とは若干異なる形となる.

Note: 信号の遅延時間が周波数によらず一定となる事は,時間軸のシフトやむだ間要素と等価であり,位相遅れが角周波数に比例するため,線形位相 (直線位相) 特性と呼ばれる.

以上の線形位相特性とは,あくまで FIR フィルタのシステム伝達関数や周波数伝達関数における特性である事に注意せよ.現実のディジタ信号処理システムでは,サンプラーの前に異名現象を低減させるローパスフィルタが設けられている場合はこれの位相特性,また 10.6 に述べた 0 次ホールドの位相特性が加わるため,これらを含めた現実のディジタル信号処理システムにおける線形位相特性の実現には高度な技法が必要となる.

13.5. 理想的ローパスフィルタのインパルス応答

式 (13‑19) を満たすインパルス応答の例として,理想的ローパスフィルタのインパルス応答を示す.このフィルタの周波数伝達関数は以下で与えられる.ここで,ωc はカットオフ周波数,T はサンプリング周期を示す.即ち「理想的」ローパスフィルタとは カットオフ周波数を超える周波数成分を全く含まないフィルタの事である.

(13-20) HdejωT= 1,ωωc0,ωc<ωπ/T

上記の周波数伝達関数と式 (12‑11) に示した離散時間フーリエ変換 (DTFT) の逆変換から,理想的ローパスフィルタのインパルス応答 hd [n] を求める.

(13-21) hdn =T2π-πTπTHdejωTejnωTdω =T2π-ωcωcejnωTdω =1πnejnωcT-e-jnωcTj2 =1πnsinnωcT,n0 hdn=T2π-ωcωcejnωTdω=ωcTπ,n=0

正弦関数は奇関数であるから,理想的ローパスフィルタのインパルス応答は hd [n] = hd [−n] となり,このフィルタの位相特性は零位相特性となる.

次に,上記のインパルス応答 hd [n] を |n| ≤ (N − 1) / 2 で打ち切り,かつ (N − 1) / 2 サンプル時間遅延させインパルス応答 h [n] を以下に示す.

(13-22) hn= sinn-N-12ωcTπn-N-12,0nN-1,nN-12ωcTπ,n=N-120,n<0,nN

上記は,式 (13‑19) を満たすため,インパルス応答を h [n] とする FIR フィルタは線形位相特性となる.

Note: 上記は N が奇数,偶数に係わらず成立する.

13.6. 線形位相 FIR フィルタの例

式 (13‑22) に示したインパルス応答に基づいた FIR フィルタのシステム伝達関数 H (z) を以下に示す.この FIR フィルタの次数は N − 1 となる.

(13-23) Hz=n=0N-1hnz-n

上記のシステム伝達関数に対応する,周波数伝達関数を以下に示す.

(13-24) HejωT=n=0N-1hne-jnωT

ωcT = π / 2 の場合の上記の周波数伝達関数のゲイン線図を図13‑5に示す.図は N = 7, 11, 19, 35, 67 の場合の例を示している.横軸の角周波数は ωT として正規化されている.また,上の図の縦軸のゲインは dB 表示,下の図の縦軸のゲインは比による表示となっている.

f:id:suzumushi0:20170611122211p:plain
図13-5: 線形位相 FIRフィルタのゲイン線図の例.

縦軸を比としたゲイン線図から明らかな通り,インパルス応答の打ち切りによって生じるギブスの現象により,ゲイン特性が波打つリップルが生じている.8.1に方形波のフーリエ級数の総和を部分和で近似した際に,ギブスの現象によって生じるリップルを示したが,上記の例の理想的ローパスフィルタのゲイン特性は,角周波数に対して方形状の特性となるため,これを近似した FIR フィルタのゲイン特性には 8.1と同様なリップルが生じる.また,縦軸を dB としたゲイン線図から明らかな通り,特に ωT > π / 2 となる周波数領域におけるゲイン特性のリップルのため,理想的ローパスフィルタと比較して,得られる減衰量に高い誤差が生じる.

このフィルタの入力の離散時間信号を x [n] = u [n] とした際の,インディシャル応答 (単位ステップ応答) y [n] を図13‑6 に示す.図から明らかな通り,このフィルタのインディシャル応答には,フィルタの次数 (レジスタ数) に比例した遅延が加わり,不連続点近傍で振動的となるリップルが生じている.

f:id:suzumushi0:20170611122212p:plain
図13-6: 線形位相 FIRフィルタのインディシャル応答の例.

13.7. 窓関数による線形位相 FIR フィルタの例

8.1 に示した通り,窓関数によってギブスの現象に起因するリップルを抑制する事が出来るため,上記の線形位相 FIR フィルタのインパルス応答にハミング窓を適用した例を示す.この例では n = (N − 1) / 2 周辺の項の重み付を高くする必要があるため,式 (8‑16) で定義した以下のハミング窓を適用する.

(13-25) wn=0.54-0.46cos2πnN-1,0nN-1

式 (13‑22) に示したインパルス応答に上記のハミング窓を適用した FIR フィルタのシステム伝達関数 H (z) を以下に示す.

(13-26) Hz=n=0N-1wnhnz-n

上記のシステム伝達関数の直流におけるゲインを 1 倍とした正規化システム伝達関数 HN (z) は以下で与えられる.

(13-27) HNz=Hzlimz1Hz=1limz1Hzn=0N-1wnhnz-n

上記の正規化システム伝達関数に対応する,周波数伝達関数を以下に示す.

(13-28) HNejωT=1limz1Hzn=0N-1wnhne-jnωT

ωcT = π / 2 の場合の上記の周波数伝達関数のゲイン線図を図13‑7に示す.横軸の角周波数は ωT として正規化されている.また,上の図の縦軸は dB 表示,下の図の縦軸は比による表示となっている.

f:id:suzumushi0:20170611122213p:plain
図13-7: 窓関数による線形位相 FIRフィルタのゲイン線図の例1.

図から明らかな通り,ハミング窓によってゲイン特性のリップルが大幅に抑制されて,減衰量が大幅に高くなっている.但し,ハミング窓による重み付の影響により,減衰傾斜は緩やかになっている.

また,このフィルタのインディシャル応答 y [n] を図13‑8に示す.図から明らかな通り,インディシャル応答には不連続点近傍で振動的となるリップルが生じている.この,リップルの幅はハミング窓が無い場合よりわずかに低くなっているが,ゲイン特性のリップルと異なり,インディシャル応答のリップルはハミング窓によって大幅には改善されない.

f:id:suzumushi0:20170611122214p:plain
図13-8: 窓関数による線形位相 FIRフィルタのインディシャル応答の例1.

以上の線形位相 FIR フィルタのゲイン線図の例は多くの教科書に示されているが,フィルタの特徴を端的に示すためにカットオフ周波数を ωcT = π / 2 ,即ちナイキスト周波数の 1 / 2 にしている様に見受けられるため,参考として ωcT = π / 10 の場合の周波数伝達関数のゲイン線図を図13‑9に示す.上の図の横軸の角周波数は ωT として正規化されている.また,下の図の横軸の角周波数は,アナログフィルタのゲイン線図の慣例に合わせ ω / ωc として正規化され,対数軸で表示されている.ここで ω / ωc = 1 は正規化されたカットオフ周波数を示す.また,この場合のインディシャル応答 y [n] を図13‑10に示す.

f:id:suzumushi0:20170611122215p:plain
図13-9: 窓関数による線形位相 FIRフィルタのゲイン線図の例2.
f:id:suzumushi0:20170611122216p:plain
図13-10: 窓関数による線形位相 FIRフィルタのインディシャル応答の例2.

本章の参考文献

A. Oppenheim, R. Schafer, 伊達訳, "ディジタル信号処理," (上), (下), コロナ社, 1978.

12. 離散フーリエ変換と離散時間フーリエ変換

本章は離散フーリエ変換 (Discrete Fourier transform: DFT) と離散時間フーリエ変換 (Discrete-time Fourier transform: DTFT) について述べる.これらは教科書に十分な説明が記載されているため,本資料では数式の証明や導出は省略し,信号の分析という観点からの解説を行う.尚,離散フーリエ変換 (DFT) と離散時間フーリエ変換 (DTFT) における種々の関係は,以下の f [n] が複素数の場合でも成立するが,その物理的な意味を見い出す事は容易ではないため,本資料では扱わない.

12.1. 離散フーリエ変換 (DFT)

離散フーリエ変換 (DFT) の導出

時間 [0, 2π] で定義された周期 2π の連続時間信号 f (t) と,9.2 に述べたサンプリング周期 T で生成される単位インパルス関数 (櫛型関数) を乗算したインパルス列を p (t) とすると,p (t), (0 ≤ t < 2π) は以下の式で記述される.尚,連続時間信号 f (t) の周期 2π 毎のサンプリング数 N = 2π / T は整数とする.

(12-1) pt=ftn=0N-1δt-nT=n=0N-1fnTδt-nT

次に,連続時間信号 f (t) とインパルス列 p (t) の関係から,離散フーリエ変換 (DFT) を導出する.式 (8‑5) に示した周期 2π の周期関数の複素フーリエ係数,及び式 (5‑7) に示した単位インパルス関数 δ (t) の定義から,p (t) の複素フーリエ係数 ck は以下で与えられる.

(12-2) ck =12π02πpte-jktdt =12π02πn=0N-1fnTδt-nTe-jktdt =12πn=0N-1fnT-ππδt-nTe-jktdt =12πn=0N-1fnTe-jknT =12πn=0N-1fnTe-j2πknN

ここで,連続時間信号 f (t) に対して,f [n] = f (nT) となる数列を離散時間信号とすると,以下の通り f [n] の離散フーリエ変換 (DFT) F [k] が得られる.尚,フーリエ変換,逆変換の定義式における係数との整合性のため F [k] = 2πck としている.

(12-3) Fk=2πck=n=0N-1fne-j2πknN

以上から明らかな通り,離散フーリエ変換 (DFT) とは本質的にインパルス列の複素フーリエ係数の別表現である

離散フーリエ変換 (DFT) の定義

f [n] = f [n + N] となる周期 N の離散時間信号 f [n] の離散フーリエ変換 (DFT) F [k],及び F [k] から f [n] を求める逆変換を以下に示す.複素指数関数は周期 2π の周期関数であるため,F [k] も周期 N の周期関数となる.以下の離散フーリエ変換 (DFT) を高速に計算するアルゴリズム高速フーリエ変換 (Fast Fourier transform: FFT) と呼ばれる

(12-4) Fk=n=0N-1fnWNkn,0kN-1 fn=1Nk=0N-1FkWN-kn,0nN-1 WN=e-j2πN

ここで,f [n] が実数,かつ N が偶数の場合,複素フーリエ級数と同様に,F [k] と F [N k], (1 ≤ k ≤ (N / 2) − 1) は共役複素数の関係となり,かつ F [0], F [N / 2] は実数となるため,逆変換の定義式は以下の通り変形できる.

(12-5) fn =1Nk=0N-1FkWN-kn=1Nk=0N-1Fkej2πknN =F0N+FN/2Nejnπ+1Nk=1N/2-1Fkej2πknN+Fke-j2πknN =F0N+FN/2Ncosnπ+1Nk=1N/2-1j2Fkej2πknN--j2Fke-j2πknNj2 =F0N+FN/2Nsinnπ+π2+k=1N/2-1Akej2πknN+φk-e-j2πknN+φkj2 Ak=2FkN,φk=Argj2FkN,1kN2-1

上記は (1) 角周波数を 2πk / N,振幅を Ak,初期位相を φk, (1 ≤ k ≤ (N / 2) − 1) とする正弦波,(2) 角周波数を π,振幅を F [N / 2] / N,初期位相を π / 2とする正弦波,及び (3) 直流成分 F [0] / N の重ね合わせの離散時間信号で f [n] を表せる事を示している.但し,周期 N / 2 はナイキスト周波数と等しいため,角周波数が π の正弦波は,振幅や初期位相によらず,振幅を F [N / 2] / N,初期位相をπ / 2とする正弦波として表現されてしまう事に注意せよ.尚,角周波数が π,初期位相が 0 となる正弦波の振幅 F [N / 2] / N は 0 と表現されるため,F [N / 2] = 0 は信号に角周波数 π の正弦波が含まれて居ない事を示している訳では無い事にも注意せよ.

Note: 高速フーリエ変換 (FFT) や逆変換,即ち離散フーリエ変換 (DFT) や逆変換は Excel のアドインで標準的に提供されており,手軽に計算する事ができる.

離散フーリエ変換 (DFT) の定理

離散フーリエ変換 (DFT) の主要な定理を以下に示す.

表12-1: 離散フーリエ変換 (DFT) の主要な定理.
線 形 性 Fafn+bgn=aFfn+bFgn=aFk+bGk
畳み込み Fm=0N-1fn-mgm=FfnFgn=FkGk
パーセバルの定理 n=0N-1fn2=1Nk=0N-1Fk2
時間軸のシフト Ffn-a=e-j2πkaNFfn=WNkaFk
周波数軸のシフト Fe-j2πnaNfn=FWNnafn=Fk+a

現実の信号分析と窓関数

高速フーリエ変換 (FFT) は,離散フーリエ変換 (DFT) を高速に計算するアルゴリズムとして信号分析において広く使われている.しかし,離散フーリエ変換 (DFT) の各種の性質が成立するのは,あくまで f [n] = f [n + N] となる周期 N の離散時間信号 f [n] に対してである.

現実の信号は周期的とは限らず,また短時間において周期的と見做せる信号,例えばモータやエンジンの回転,であっても,負荷の変動による駆動系の捩れ等により,周期の遅れや進みが生じるため,その周期に合わせて離散時間信号をサンプリングする事は容易では無い.そもそも信号にどの様な周波数成分が含まれているか判らないから分析するのであって,最初から信号が厳密に周期的と判っているのであれば,分析する必要も無いのである

従って,現実の高速フーリエ変換 (FFT) による信号分析とは,離散時間信号を N 回サンプリングし,それを周期 N の離散時間信号であると仮定して離散フーリエ変換 (DFT) を計算する事となる.この場合,離散時間信号に非周期信号や周期が N とはならない周期信号や含まれている可能性は十分にあり,これらは離散フーリエ変換 (DFT) の計算結果に角周波数や振幅,初期位相が変動する信号やノイズとして現れる.

これらの分析結果の変動やノイズを低減する技法として,8.1 に述べた窓関数が広く用いられている.代表的な窓関数として,式 (8-16) に示したハミング窓を以下に再掲する.

(12-6) wn=0.54-0.46cos2πnN-1,0nN-1

上記のハミング窓を離散時間信号 f [n] と乗算し w [n] f [n] とすると,ハミング窓の n = N / 2 を中心とした項の重みと比較して,n = 0 及び n = N − 1 周辺の項の重みは非常に低いため,w [0] f [0] と w [N − 1] f [N − 1] の周辺の差が減少し,w [n] f [n] を連続的,周期的な信号と見做せるようになる.このため,w [n] f [n] の離散フーリエ変換 (DFT) を計算すれば分析結果の変動やノイズを低減する事ができる.その反面,得られる分析結果の精度は低下する.

12.2. 離散時間フーリエ変換 (DTFT)

離散時間フーリエ変換 (DTFT) の導出

時間 [−∞, ∞] で定義された連続時間信号 f (t) と,9.2 に述べたサンプリング周期 T で生成される単位インパルス関数 (櫛型関数) を乗算したインパルス列を p (t) とすると,p (t) は以下の式で記述される.

(12-7) pt=ftn=-δt-nT=n=-fnTδt-nT

次に,連続時間信号 f (t) とインパルス列 p (t) の関係から,離散時間フーリエ変換 (DTFT) を導出する.式 (8‑19) に示したフーリエ変換の定義,及び式 (5‑5) に示した単位インパルス関数 δ (t) の定義から,p (t) のフーリエ変換 P (ω) は以下で与えられる.

(12-8) Pω =-pte-jωtdt =-n=-fnTδt-nTe-jωtdt =n=-fnT-δt-nTe-jωtdt =n=-fnTe-jnωT

ここで,連続時間信号 f (t) に対して,f [n] = f (nT) となる数列を離散時間信号とすると,以下の通り離散時間フーリエ変換 (DTFT) が得られる.

(12-9) Pω=n=-fne-jnωT

以上から明らかな通り,離散時間フーリエ変換 (DTFT) とはインパルス列のフーリエ変換の別表現である

離散時間フーリエ変換 (DTFT) の定義

サンプリング周期を T とする離散時間信号 f [n] の離散時間フーリエ変換 (DTFT) F (ω) を以下に示す.複素指数関数は周期 2π の周期関数であるため,F (ω) は周期 2π / T の周期関数となる.

(12-10) Fω=n=-fne-jnωT

F (ω) から f [n] を求める逆変換を以下に示す.

(12-11) fn=T2π-π/Tπ/TFωejnωTdω

ここで,f [n] が実数の場合は,フーリエ逆変換と同様に,F (ω) と F (−ω), (0 < ω ≤ (π / T)) は共役複素数の関係となるため,逆変換の定義式 (12‑11) は以下の通り変形できる.

(12-12) fn=T2π-πTπTFωejnωTdω fn-TF02π ~T2π0+πTFωejnωT+FωejnωTdω =Tπ0+πTjFωejnωT--jFωejnωTj2dω =0+πTAωejnωT+φω-e-jnωT+φωj2dω Aω=TFωπ,φω=ArgjTFωπ,0<ωπT

上記は角周波数をω,振幅を A (ω),初期位相を φ (ω), (0 < ω ≤ (π / T)) とする正弦波と,直流成分 T F (0) / (2π) の重ね合わせの離散時間信号で f [n] を表せる事を示している.

Note: 離散時間フーリエ変換 (DTFT) はインパルス列のフーリエ変換の別表現に過ぎないため,式 (8‑21) に示したフーリエ逆変換によって,F (ω) からインパルス列を得る事もできる.

フーリエ変換離と散時間フーリエ変換 (DTFT) の関係

連続時間信号 f (t) のフーリエ変換Fa (ω),連続時間信号 f (t) を周期 T でサンプリングした離散時間信号を f [n],この離散時間フーリエ変換 (DTFT) を F (ω) とすると,F (ω) と Fa (ω) には以下の関係がある.(下記の導出は煩雑なため後述する).

(12-13) Fω=1Tn=-Faω-2πnT

上記の関係は 9.1 に述べたサンプリングにおける異名現象 (エイリアシング) の発生を示している.また,上記は F (ω) が周期 2π / T の周期関数となる事,及び 10.2 に述べたインパルス列による階段関数の信号の近似により,離散時間フーリエ変換 (DTFT) はフーリエ変換と比較して1 / T 倍のゲインが生じる事も示している.

フーリエ級数展開と離散時間フーリエ変換 (DTFT) の関係

以上述べた通り,散時間フーリエ変換 (DTFT) は,離散時間の信号を連続かつ周期的な角周波数の正弦波と直流成分の重ね合わせで表す.一方,8.1 に述べたフーリエ係数は,周期的な連続時間の信号を離散的な角周波数の正弦波と直流成分の重ね合わせで表す.このため,式 (8‑12) に示した周期 2π / ω の周期関数の複素フーリエ係数と級数は,式 (12‑10) 及び (12‑11) に示した離散時間フーリエ変換 (DTFT),逆変換と,時間と角周波数を入れ替えた対称的な関係にある.

具体的には,複素フーリエ級数 f (t) を離散時間フーリエ変換 (DTFT) F (ω) に,複素フーリエ係数 cn を逆変換 f [n] に,複素フーリエ級数,係数の ω, t を各々離散フーリエ変換,逆変換の T, ω に置き換えると,これらは実質的に同じ式となる.このため,フーリエ級数展開における各種の性質,例えばギブスの現象は,時間と角周波数を入れ替えた関係で,離散時間フーリエ変換 (DTFT) においても発生し,これによって生じる周波数特性のリップルを抑制する技法として窓関数が使用できる点も同じとなる

離散時間フーリエ変換 (DTFT) の定理

離散時間フーリエ変換 (DTFT) はフーリエ変換の別表現であるから,これと同様な定理が成立する.離散時間フーリエ変換 (DTFT) の主要な定理を以下に示す.

表12-2: 離散時間フーリエ変換 (DTFT) の主要な定理.
線 形 性 Fafn+bgn=aFfn+bFgn=aFω+bGω
畳み込み Fk=-fn-kgk=FfnFgn=FωGω
パーセバルの定理 n=-fn2=12π-π/Tπ/TFω2dω
時間軸のシフト Ffn-a=e-jaωTFfn=e-jaωTFω
周波数軸のシフト Fe-jatfn=Fω+a

参考: フーリエ変換と離散時間フーリエ変換 (DTFT) の関係の導出

式 (12-13) に示したフーリエ変換と離散時間フーリエ変換 (DTFT) の関係式の導出を以下に示す.この導出は煩雑であるため,ステップに別けて示す.

フーリエ変換における周波数関数の畳み込み

表8-1に示した,フーリエ変換の畳み込みの定理を以下に示す.

(12-14) F-ft-τgτdτ=FωGω

上記は,時間関数の畳み込みであるが,フーリエ変換の対称性より,以下の通り周波数関数に対しても畳み込みが成立する.

(12-15) F12π-Fω-τGτdτ=ftgt

櫛型関数のフーリエ級数

9.2 のサンプラーの論理モデルにおいて述べた,サンプリング周期 T で生成される単位インパルス関数 (櫛型関数) のフーリエ級数を示す.櫛型関数は周期 T の周期関数であるから,これのフーリエ係数 Cn は以下で与えられる.以下の展開において,積分範囲において 0 以外となる単位インパルス関数は 1 つしか無い事を利用している.

(12-16) cn =1T-T/2T/2n=-δt-nTe-j2πntTdt =1Tn=--T/2T/2δt-nTe-j2πntTdt =1Tn=--T/2T/2δte-j2πntTdt =1Tn=-e0=1T

従って,櫛型関数のフーリエ級数は以下で与えられる.

(12-17) n=-δt-nT=1Tn=-ej2πntT

櫛型関数のフーリエ変換

離散時間フーリエ変換 (DTFT) はインパルス列のフーリエ変換の別表現であるから,櫛型関数のフーリエ変換は,以下に示す通りf [n] = 1 とした離散時間フーリエ変換 (DTFT) の定義式から与えられる.

(12-18) Fn=-δt-nT=n=-e-jnωT

しかし,上に述べた櫛型関数のフーリエ級数,及び式 (8-28) に示した単位インパルス関数の積分表示から,櫛型関数のフーリエ変換は,以下の別形式で表す事もできる.

(12-19) Fn=-δt-nT =-n=-δt-nTe-jωtdt =1T-n=-ej2πntTe-jωtdt =1Tn=--ej2πntTe-jωtdt =1Tn=--e-jω-2πnTtdt =2πTn=-δω-2πnT

フーリエ変換と離散時間フーリエ変換 (DTFT) の関係の導出

以上から式 (12‑13) を導出する.連続時間信号 f (t) を周期 T でサンプリングしたインパルス列 p (t) は以下で与えらえる.

(12-20) pt=ftn=-δt-nT

上記は連続時間信号と櫛型関数の積であるから,インパルス列 p (t) のフーリエ変換,即ち離散時間信号 f [n] の離散時間フーリエ変換 (DTFT) F (ω) は,以下の通り,連続時間信号のフーリエ変換 Fa (ω) と式 (12‑19) に示した櫛型関数のフーリエ変換の畳み込みで与えられる.

(12-21) Fω =12π-Faτ2πTn=-δω-2πnT-τdτ =1Tn=--Faτδω-2πnT-τdτ =1Tn=-Faω-2πnT

本章の参考文献

A. Oppenheim, R. Schafer, 伊達訳, "ディジタル信号処理," (上), (下), コロナ社, 1978.
E. Oran Brigham, 宮川, 今井訳, "高速フーリエ変換," 科学技術出版社, 1985.

11. IIR フィルタの概要

本章では,11.1 において IIR (Infinite Impulse Responce) フィルタ (無限インパルス応答フィルタ,巡回型フィルタ) の構造,11.2 においてインパルス不変設計法,11.3 において双一次変換による設計法を復習する.これらは教科書に十分な説明が記載されているため,本資料は留意すべき点を Note に示す.

11.1. IIR フィルタの構造

システム伝達関数 H (z) が以下の式で与えられるディジタルフィルタの構造を示す.ここで,X (z), Y(z) は各々ディジタルフィルタの入出力となる離散時間信号 x [n], y [n], (n ≥ 0) の Z 変換を示す.

(11-1) Hz=YzXz=r=0Mbrz-r1-k=1Nakz-k

上記のディジタルフィルタに入出力される離散時間信号 x [n], y [n] の関係は,以下の差分方程式で記述される.

(11-2) yn=k=1Nakyn-k+r=0Mbrxn-r

IIRフィルタとは,上記の差分方程式を,フィルタの実現アルゴリズムと見做したものであり,出力の離散時間信号 y [n] を,過去 N 回の出力 y [n − 1], ..., y [n N] と,現在及び過去 M 回の入力 x [n], ..., x [n M] から求める.

基本的な IIR フィルタの構造を図11‑1 に示す.図において,四角は入力信号に 1 サンプリング周期の遅延を加えて出力するレジスタ (Z 変換において1 サンプリング周期の時間シフトは z−1 となる事からこの様に表現される),三角は定数係数の乗算器,丸は加算器を示す.レジスタの初期値は全て 0 となる.また,図を簡略化するために M = N としている.N, M は各々システム伝達関数の分子の次数,分母の次数と呼ばれる.

図11-1 の左の図は,最も基本的な形の IIR フィルタで,上記の差分方程式を式の通りに計算する.また,現在及び過去の入力の計算と,過去の出力の計算は入れ替える事が出来るため,これを右の図の様に変形させて,レジスタ数を削減する場合が多い.これらは直接形の IIR フィルタと呼ばれる.

f:id:suzumushi0:20170611122159p:plain
図11-1: IIR フィルタの構造 (直接形).

また,システム伝達関数 H (z) は,以下の様に因数分解する事ができる.ここで,⌊x⌋ は x を超えない最大の整数を示す.

(11-3) Hz=r=0Mbrz-r1-k=1Nakz-k=Ak=1N+1/21+β1kz-1+β2kz-21-α1kz-1-α2kz-2,MN

因数分解された各項を直接形の IIR フィルタで構成し,以下の図に示す様に,これらを縦続接続する事によって,システム伝達関数で記述された特性のフィルタを得る事もできる.これは,縦続形の IIR フィルタと呼ばれる.

f:id:suzumushi0:20170611122200p:plain
図11-2: IIR フィルタの構造 (縦続形).

11.2. インパルス不変設計法

インパルス不変設計法とは,アナログフィルタの伝達関数 G (s) から,ディジタルフィルタのシステム伝達関数 H (z) を求める設計法であり,以下に示す様に,ディジタルフィルタのインパルス応答の離散時間信号 h [n], (n ≥ 0) は,アナログフィルタのインパルス応答 g (t), (t ≥ 0) のサンプリングとなる.

(11-4) hn=gnT,n0

Note: 上記は教科書的なインパルス不変設計法の定義であり,10 章に述べたゲインや誤差の補正は考慮されていない.また,9.10 に述べたローパスフィルタの設計例は,この教科書的な定義に基づいている.

システム伝達関数

以下の伝達関数で示される,時定数が τ,カットオフ周波数が ωc = 1 / τ のローパスフィルタを,インパルス不変設計法に基づいてディジタル化した例を示す.

(11-5) Gs=1τs+1=ωcs+ωc

表5‑1 より,上記の伝達関数で表されるローパスフィルタのインパルス応答 g (t) は以下の通りとなる.

(11-6) gt=1τe-tτ,t0

また,表9‑3 より,上記の伝達関数に対応するシステム伝達関数を以下に示す.サンプリング周期は T とする.

(11-7) Hz=ωc1-e-ωcTz-1

表9‑1 より,上記のシステム伝達関数で表されるローパルスフィルタのインパルス応答 h [n] は以下の通りとなる.このインパルス応答は,明らかに,式 (11‑6) に示した連続時間のインパルス応答 g (t) を離散時間化したものとなる.

(11-8) hn=1τe-nTτ,n0

正規化システム伝達関数と周波数伝達関数

上記のシステム伝達関数のゲインや誤差を補正した,正規化システム伝達関数は以下で与えられる.

(11-9) HNz=lims0Gslimz1HzHz=1-e-ωcT1-e-ωcTz-1

上記の正規化システム伝達関数に対応する,周波数伝達関数を以下に示す.

(11-10) HNejωT=1-e-ωcT1-e-ωcTe-jωT

差分方程式

正規化システム伝達関数は,ディジタルフィルタの入出力となる離散時間信号 x [n], y [n] の Z 変換の比であるから,これを以下の様に逆 Z 変換する事によって,離散時間信号 x [n], y [n] の関係を記述する差分方程式が得られる.

(11-11) HNz=YzXz=1-e-ωcT1-e-ωcTz-1 1-e-ωcTz-1Yz=1-e-ωcTXz yn-e-ωcTyn-1=1-e-ωcTxn yn=e-ωcTyn-1+1-e-ωcTxn

ボード線図

上記の周波数伝達関数のボード線図を図11‑3に示す.横軸の角周波数は ω / ωc として正規化されている.図中の青線は正規化されたナイキスト周波数 fn を示しており,この例ではカットオフ周波数 ωc の 100 倍の角周波数としている.また,図中の赤線はアナログフィルタの特性を示す.

f:id:suzumushi0:20170611122201p:plain
図11-3: インパルス不変設計法によるローパスフィルタのボード線図の例.

図から明らかな通り,ナイキスト周波数がカットオフ周波数と比較して十分に高い場合は,この設計法によるローパスフィルタのゲイン特性は,ナイキスト周波数近傍で減衰量が低下している以外は,アナログフィルタの特性と高い精度で一致する.一方,ω / ωc ≪ 1 の場合の位相特性はアナログフィルタの特性と高い精度一致するが,ω / ωc ≫ 1 の場合はかなりの相違がみられる.

ナイキスト周波数をカットオフ周波数と比較して十分に高く設計する事は,経済的では無い場合が多いため,次にナイキスト周波数をカットオフ周波数 ωc の 20 倍の角周波数とした場合のボード線図を示す.

f:id:suzumushi0:20170611122202p:plain
図11-3: インパルス不変設計法によるローパスフィルタのボード線図の例.

この場合のナイキスト周波数以下のゲイン特性は,前の例と同様に,ナイキスト周波数近傍で減衰量が低下している以外は,アナログフィルタの特性と高い精度で一致する.

また,9.1 に述べた通り,サンプリング周波数 fs = 1 / T の離散時間信号では,異名現象により,ナイキスト周波数 fn = fs / 2 以下の周波数成分 f0 と,ナイキスト周波数を超える周波数成分 N fs ± f0, (N = 1, 2, 3, ...) を区別する事はできないため,図に示した通り,ナイキスト周波数を超えるゲイン特性には折り返しが生じる.

ナイキスト周波数近傍で減衰量が低下する原因も異名現象に起因する.バタワース,チェビシェフ,ベッセル,楕円等の既存のアナログローパスフィルタのゲイン特性は,或る特定の角周波数,若しくは角周波数 ∞ においてゲインが 0 倍となる特性であり,或る角周波数以上においてゲインが 0 倍となる特性では無いため,これらのフィルタではナイキスト周波数以上において,必ず或る程度のゲインが生じる.このため,下図に示す様に,異名現象によるゲイン特性の折り返し (重複) により,特にナイキスト周波数付近の減衰量が低下する.

f:id:suzumushi0:20170611122203p:plain
図11-5: インパルス不変設計法によるローパスフィルタのゲイン特性における異名現象.

インパルス不変設計法によるローパスフィルタでは,この異名現象による減衰量の低下は大なり小なり生じるため,設計時にそれが許容範囲内にある事を確認する必要がある.

11.3. 双一次変換による設計法

双一次変換による設計法とは,アナログフィルタの伝達関数 G (s) から,ディジタルフィルタのシステム伝達関数 H (z) を求める設計法であり,インパル不変設計法における減衰量の低下を抑える事ができる.

11.2 に述べた異名現象によるゲイン特性の折り返しを避けるため,双一次変換による設計法では,アナログフィルタの周波数 0 から ∞ までの特性が,周波数 0 からナイキスト周波数 fn までの特性となる様に周波数軸を変換したフィルタを考える.この様な周波数軸の写像は無数に存在するが,双一次変換では,以下に示すアナログフィルタの角周波数ωa と,周波数軸が変換されたフィルタの角周波数 ωの関係が用いられる.ここで T はサンプリング周期である.

(11-12) ωa=tanωT2

上記の関係ではアナログフィルタの角周波数 ωa = 0 はω = 0 に,ωa = ∞ はω = π / T = π fs = 2 π fn 即ちナイキスト周波数写像される.この関係に基づいて周波数軸が変換されたフィルタをディジタル化すると,以下の図に示す様に,異名現象は生じるが,ゲイン特性の折り返し (重複) は生じないため,インパルス不変設計法において見られた減衰量の低下は生じない.但し,その代償として,上記の角周波数の対応関係の歪により,特に角周波数が高い領域において,アナログフィルタと周波数特性が異なるディジタルフィルタが得られる.

f:id:suzumushi0:20170611122204p:plain
図11-6: 双一次変換による設計法の原理.

上記の角周波数の対応関係は,以下の様に変形できる.

(11-13) jωa=jtanωT2=jsinωT2cosωT2=jejωT2-e-jωT22jejωT2+e-jωT22=ejωT-1ejωT+1

ディジタル信号処理では,初期条件を 0 とする事が暗黙の前提であるため,上記のasa に,s に各々置き換える事によって,以下に示す通り,アナログフィルタの伝達関数 G (sa) と,周波数軸を変換したフィルタの伝達関数 G (s) の対応関係が得られる.

(11-14) sa=esT-1esT+1

Z 変換とはラプラス変換sz = eTs に置き換えたものであるから,双一次変換では,アナログフィルタの伝達関数G (s) の sa を下記の式の通り z に置き換え,ディジタルフィルタのシステム伝達関数 H (z) とする.

(11-15) sa=esT-1esT+1=z-1z+1=1-z-11+z-1

ローパスフィルタのシステム伝達関数と周波数伝達関数

以下の伝達関数で示される,時定数が τ,カットオフ周波数が ωc = 1 / τ のローパスフィルタを,双一次変換による設計法に基づいてディジタル化した例を示す.

(11-16) Gs=1τs+1=ωcs+ωc

この設計法によるアナログフィルタとディジタルフィルタの角周波数の対応関係には歪が生じるため,予め上記のアナログフィルタの時定数を補正する必要がある.式 (11‑12) に基づいて補正したカットオフ周波数周波数ωa を以下に示す.サンプリング周期は T とする.

(11-17) ωa=tanωcT2

上記の伝達関数に対応するシステム伝達関数を以下に示す.

(11-18) Hz=ωas+ωa=ωa1-z-11+z-1+ωa=ωa1+z-11+ωa-1-ωaz-1

また,上記のシステム伝達関数に対応する,周波数伝達関数を以下に示す.

(11-19) HejωT=ωa1+e-jωT1+ωa-1-ωae-jωT

Note: 双一次変換は本質的に周波数軸の変換であるため,上記のシステム伝達関数のゲインを補正する必要は無い.(上記の周波数伝達関数ω → 0,即ち直流におけるゲインは 1 倍であり,アナログフィルタのゲインと一致している).

ローパスフィルタの差分方程式

インパルス不変設計法と同様に,上記のシステム伝達関数で記述されるローパスフィルタの差分方程式は以下の通りとなる.

(11-20) Hz=YzXz=ωa1+z-11+ωa-1-ωaz-1 1+ωa-1-ωaz-1Yz=ωa1+z-1Xz 1+ωayn-1-ωayn-1=ωaxn+xn-1 yn=1-ωa1+ωayn-1+ωa1+ωaxn+xn-1

ローパスフィルタのボード線図

上記の周波数伝達関数のボード線図を図11‑7に示す.横軸の角周波数は ω / ωc として正規化されている.図中の青線は正規化されたナイキスト周波数 fn を示しており,この例ではカットオフ周波数 ωc の 20 倍の角周波数としている.また,図中の赤線はアナログフィルタの特性を示す.

f:id:suzumushi0:20170611122205p:plain
図11-7: 双一次変換による設計法におけるローパスフィルタのボード線図の例.

図から明らかな通り,この設計法によるローパスフィルタのゲイン特性は,アナログフィルタとディジタルフィルタの角周波数の対応関係における歪のため,角周波数が高い領域において減衰量が急上昇している以外は,アナログフィルタの特性と高い精度で一致する.また,異名現象によりナイキスト周波数を超えるゲイン特性には折り返しが生じているが,インパルス不変設計法に見られたナイキスト周波数近傍における減衰量の低下は生じていない.

ハイパスフィルタのシステム伝達関数と周波数伝達関数

以下の伝達関数で示される,時定数が τ,カットオフ周波数が ωc = 1 / τ のハイパスフィルタを,双一次変換による設計法に基づいてディジタル化した例を示す.

(11-21) Gs=τsτs+1=1-1τs+1=1-ωcs+ωc

ローパスフィルタと同様に,予め上記のアナログフィルタの時定数を補正する必要がある.式 (11‑12) に基づいて補正したカットオフ周波数 ωa を以下に示す.サンプリング周期は T とする.

(11-22) ωa=tanωcT2

上記の伝達関数に対応するシステム伝達関数を以下に示す.

(11-23) Hz=1-ωas+ωa=1-ωa1-z-11+z-1+ωa=1-ωa1+z-11+ωa-1-ωaz-1

また,上記のシステム伝達関数に対応する,周波数伝達関数を以下に示す.

(11-24) HejωT=1-ωa1+e-jωT1+ωa-1-ωae-jωT

ハイパスフィルタの差分方程式

上記のシステム伝達関数で記述されるハイパスフィルタの差分方程式は以下の通りとなる.

(11-25) Hz=YzXz=1-ωa1+z-11+ωa-1-ωaz-1 1+ωa-1-ωaz-1Yz=1-z-1Xz 1+ωayn-1-ωayn-1=xn-xn-1 yn=1-ωa1+ωayn-1+11+ωaxn-xn-1

ハイパスフィルタのボード線図

上記の周波数伝達関数のボード線図を図11‑8に示す.横軸の角周波数は ω / ωc として正規化されている.図中の青線は正規化されたナイキスト周波数 fn を示しており,この例ではカットオフ周波数 ωc の 5 倍の角周波数としている.また,図中の赤線はアナログフィルタの特性を示す.

f:id:suzumushi0:20170611122206p:plain
図11-8: 双一次変換による設計法におけるハイパスフィルタのボード線図の例.

図から明らかな通り,この設計法によるハイパスフィルタのナイキスト周波数以下のゲイン特性は,ナイキスト周波数 fn がカットオフ周波数 ωc の僅か 5 倍にも係わらず,アナログフィルタの特性と高い精度で一致する.ローパスフィルタと同様に,角周波数が高い領域において,アナログフィルタとディジタルフィルタの角周波数の対応関係に歪は生じているが,ローパスフィルタでは角周波数が高い領域のゲインは一定であるため,結果として影響は無視できる.一方,位相特性には,角周波数の対応関係の歪による影響が明確に現れている.

参考: インパルス不変設計法の適用範囲

一部の教科書には,ハイパスフィルタ,及びノッチフィルタ (帯域消去フィルタ) は,ナイキスト周波数近傍での帯域制限が困難であり,異名現象によるゲイン特性への影響度が高いため,インパルス不変設計法は,これらのフィルタの設計には適さないとされている.本資料に述べた通り,インパルス不変設計法では,ローパスフィルタにおいても異名現象による減衰量の低下は大なり小なり生じるため,それを許容できるかは単なる程度問題でしか無い.このため,インパルス不変設計法によるハイパスフィルタ,及びノッチフィルタの特性を不適切とする根拠となる設計例を探したが見当たらなかった.

そもそも,異名現象より遥か以前の問題として,インパルス不変設計法によるハイパスフィルタ及びノッチフィルタの設計自体が困難である.ハイパスフィルタの伝達関数の例を以下に示す.

(11-26) Gs=τsτs+1=1-1τs+1

上記のインパルス応答 g (t) を以下に示す.

(11-27) gt=δt-1τe-tτ,t0

上記から明らかな通り,ハイパスフィルタ及びノッチフィルタにおける伝達関数の分子と分母の次数は等しいため,5.4 に述べた通り,これらのインパルス応答には必ず単位インパルス関数が現れる.これを 9.2 に述べた櫛型関数によってサンプリングした場合,インパルス列に単位インパルス関数の 2 乗が現れ,離散時間信号 h [0] は発散する

従って,ディジタル化されたハイパスフィルタ,及びノッチフィルタのインパルス応答の離散時間信号を,厳密にアナログフィルタのインパルス応答のサンプリングとなる様に設計する事は不可能である.但し,単位インパルス関数の離散時間信号を 1 や 1 / T 等とするご都合主義的な解釈も考えられ,それに基づいてハイパスフィルタらしきものはできるが,その特性を正確に議論する事は困難となる.

上記の様にインパルス応答に単位インパルス関数が現れる場合は,インディシャル応答 (単位ステップ応答) 不変設計法により,適切な特性のフィルタを設計できるが,双一次変換による設計法と比較して特に優位性のある結果は得られないため,詳細は割愛する.

尚,インパルス不変設計法に基づいてローパスフィルタを設計し,これをハイパスフィルタに変換した場合は,変換元のローパスフィルタに生じている異名現象による特性も含めて,変換先のハイパスフィルタの特性に変換されるため,不適切とするほどの異名現象による影響は生じない.

本章の参考文献

A. Oppenheim, R. Schafer, 伊達訳, "ディジタル信号処理," (上), (下), コロナ社, 1978.

10. ディジタル信号処理の論理モデルと実装モデル

本章は,先ず,10.1において,ディジタル信号処理の理論モデルと実装モデルを述べ,システム伝達関数は,物理的に実存しない想像上の概念であるインパルス列の関係を記述している事を示す.次に,10.2 において,システム伝達関数を実装モデルに適用するために,論理モデルのインパルス列と実装モデルの階段関数の信号との関係に基づいたゲインの補正方法を述べ,10.3において,この方法では誤差が生じる事を示す.このため10.4 及び 10.5において,アナログフィルタとディジタルフィルタにおけるインパルス応答の関係に基づいた,この誤差の補正方法を明らかにする.最後に10.6において 0 次ホールドの伝達関数を示し,ディジタルフィルタの設計には,これも考慮する必要がある事を述べる.

10.1. ディジタル信号処理の論理モデルと実装モデルの関係

ディジタル信号処理の論理モデルと実装モデルを図10-1 に示す.上のブロックダイヤグラムは論理モデルを,下のブロックダイヤグラムは実装モデルを表す.

f:id:suzumushi0:20170611122152p:plain
図10-1: ディジタル信号処理の論理モデルと実装モデル.

ディジタル信号処理の論理モデル

9.2 に述べた通り,ディジタル信号処理の論理モデルでは,サンプラーにおいて,連続時間信号の入力 x (t) はサンプリング周期 T で生成される単位インパルス関数 (櫛型関数) と乗算され,インパルス列 p (t) に変換される.インパルス列 p (t) における各単位インパルス関数の係数は離散時間信号 x [n] となり,ディジタル信号処理の入力となる.

そして,ディジタル信号処理の出力の離散時間信号 y [n] は,再びインパルス列 q (t) と見做され 0 次ホールドに入力される.0 次ホールドはインパルス列 q (t) を階段関数の信号 y (t) に変換して出力する.

ディジタル信号処理の実装モデル

一方,ディジタル信号処理の実装モデルでは,連続時間信号の入力 x (t) は最初に SH (Sample and hold) に入力される.これは,次段の A/D 変換器 (Analog-to-digital converter) が,アナログ入力からディジタル出力への変換に要する時間の間,A/D 変換器の入力を安定化させる目的で設けられている.SH は,次段の A/D 変換器がアナログ入力をディジタル出力に変換する際に,連続時間信号の入力 x (t) をごく短時間保持して A/D 変換器に出力する.尚,これ以外のタイミングにおいて SH は x (t) をそのまま出力する.

A/D 変換器は,サンプリング周期 T 毎に,アナログ入力をディジタル出力 (量子化された離散時間信号 x [n]) に変換する.ディジタル信号処理部分は,論理回路DSP,コンピュータ等で構成される.D/A 変換器 (Digital-to-analog converter) は,サンプリング周期 T 毎に,ディジタル入力 (量子化された離散時間信号 y [n]) をアナログ出力 (階段関数の信号 y (t)) に変換する.

出力側の SH は,前段の D/A 変換器のディジタル入力が変化する際に生じるノイズを避ける目的で設けられている.多くの場合,D/A 変換器ではディジタル入力の MSB が変化する際に大きなノイズが生じる.このため,SH は前段の D/A 変換器のディジタル入力が変化する際に,出力信号 y (t) をごく短時間保持する.尚,これは簡易な用途では省略される事が多い.

論理モデルと実装モデルの関係

何れのモデルも,連続時間信号の周期的な瞬時値を離散時間信号としている点は同じである.しかし,これらのモデルにおける離散時間信号の意味は全く異なる.論理モデルにおける離散時間信号とは,インパルス列における各単位インパルス関数の係数の事である.一方,実装モデルにおける離散時間信号は,図10-2 に示す様に D/A 変換器の出力が階段関数となる事から,階段関数の信号における各階段の高さを示す係数を意味する.

f:id:suzumushi0:20170611122153p:plain
図10-2: 階段関数の信号.

論理モデルにおけるサンプラーは,連続時間信号の入力 x (t) を,インパルス列 p (t) に変換する.一方,実装モデルにおける SH と A/D 変換器は,連続時間信号の入力 x (t) を離散時間信号 x [n] に変換する.また,論理モデルにおける 0 次ホールドは,インパルス列 q (t) を階段関数の信号 y (t) に変換する.一方,実装モデルにおける D/A 変換器と SH は,離散時間信号 y [n] を階段関数の信号 y (t) に変換する.

Note: SH は,サンプリング周期毎に,入力信号をごく短時間保持して出力する機能であり,論理モデルにおけるサンプラーとは全く異なる.

以上から明らかな通り,ディジタル信号処理の実装モデルにはインパルス列は存在しない.単位インパルス関数は物理的には実存しない理論上,想像上の関数であるのと同様に,インパルス列も理論上,想像上の概念である.一方,9.9 に述べた通り,システム伝達関数はシステムに入出力されるインパルス列のラプラス変換の比,即ち Z 変換の比を記述する.即ち,ディジタル信号処理の論理モデルは,実存しない想像上の信号の関係を論じているのである.従って,システム伝達関数がシステムに入出力される離散時間信号の比を記述していると誤解して,ディジタルフィルタを実装しても期待した特性のフィルタは得られないのである.

Note: 非常に不思議な事に多くのディジタル信号処理の教科書では,本資料の9.2サンプル値の表現,及び9.3 Z 変換の導出の部分が省略されている.9.4 に示した通り Z 変換は離散時間信号を元に定義されているため,これらの知識が無いと,Z 変換は離散時間信号の母関数の一種であると理解してしまい,システム伝達関数はシステムに入出力される離散時間信号の比を記述していると誤解するのである.

10.2. インパルス列と階段関数の信号の関係

システム伝達関数をディジタル信号処理の実装モデルに適用するために,同一の離散時間信号で表示される,論理モデルのインパルス列と実装モデルの階段関数の信号との関係を示す.階段関数の信号とは連続時間信号の長方形近似であるから,論理モデルのインパルス列を実装モデルの階段関数の信号の近似と考える.ここで,係数が x [n] となる,階段関数の中の一つの長方形と,インパルス列中の一つのインパルスを図10‑3 に示す.サンプリング周期は T とする.

f:id:suzumushi0:20170611122154p:plain
図10-3: 長方形とインパルス.

この近似の精度は長方形とインパルスの面積比から与えられる.単位インパルス関数の定義より,上記のインパルスの面積は以下となる.

(10-1) xn-δt-nTdt=xn

一方,上記の長方形の面積は T x [n] であるから,階段関数の信号を,これと同一の離散時間信号で表示されるインパルス列で近似すると 1 / T 倍の「近似」となる.即ち,ある離散時間信号で表示された階段関数の信号を,同じ離散時間信号で表示されたインパルス列に変換すると,1 / T 倍のゲインが生じる.同様に,ある離散時間信号で表示されたインパルス列を,同じ離散時間信号で表示された階段関数の信号に変換すると T 倍のゲインが生じる.

従って,図10‑1 に示したディジタル信号処理の実装モデルにおける離散時間信号 x [n] は階段関数の信号の係数であるから,システム伝達関数を実装モデルに適用する場合は,これを T 倍する必要がある.同様に,ディジタル信号処理の論理モデルにおける離散時間信号 y [n] はインパルス列 q (t) の係数であるから,0 次ホールドの伝達関数を実装モデルに適用する場合は,これを1 / T 倍する必要がある.

これが,9.10 に示したローパスフィルタの例において,サンプリング周波数を 1 [MHz] = 106 [Hz] としたため,ゲインの最大値が 120 dB 以上となった理由である.また,ディジタル信号処理の教科書における離散時間信号やシステム伝達関数等の式の係数が微妙に異なる理由でもある.

10.3. ゲインを補正したローパスフィルタの例

9.10 に示したローパスフィルタのシステム伝達関数は,上に述べた論理モデルのインパルス列と実装モデルの階段関数の信号とのゲインの関係に基づいて,以下の通り補正される.

(10-2) THz=ωcT1-e-ωcTz-1

また,上記のゲインが補正されたシステム伝達関数に対応する,周波数伝達関数を以下に示す.

(10-3) THejωT=ωcT1-e-ωcTe-jωT

上記の周波数伝達関数のゲイン線図を図10‑4に示す.横軸の角周波数は ω / ωc として正規化されている.また図中の青線は正規化されたナイキスト周波数 fn を示しており,この例ではカットオフ周波数 ωc の 20 倍の角周波数としている.

f:id:suzumushi0:20170611122155p:plain
図10-4: ゲインを補正ローパスフィルタのゲイン線図の例.

上記のゲイン線図と図6-8 に示した一次遅れ要素 (ローパスフィルタ) のゲイン線図を比較すると,ナイキスト周波数以下の特性は一見似ている様に見えるが,ゲインの最大値が 0 dB 以上 (この例では約 0.7 dB,約 1.08 倍) となっており特性が微妙に異なるローパスフィルタとなっている.

ここで,ゲインが補正された周波数伝達関数における,直流のゲインを以下に示す.

(10-4) limω0THejωT=limω0ωcT1-e-ωcTe-jωT=ωcT1-e-ωcT

上記の式の指数関数をテーラー展開すれば,ゲインが約 1 倍となる事は判るが,明らかに厳密に 1 倍とはならない.この誤差の主な要因は,実装モデルでは階段関数の信号によって,連続時間信号を近似している点にある.(根拠は後述する).階段関数によって近似された連続時間信号の例を図10-5 に示す.

f:id:suzumushi0:20170611122156p:plain
図10-5: 階段関数による連続時間信号の近似.

図に示した様な下に凸な関数を階段関数で近似した場合,元の関数より過大な関数となる.同様に上に凸な関数を階段関数で近似した場合,元の関数より過小な関数となる.

10.4. 階段関数の信号による誤差の補正

そこで,アナログフィルタに基づいて設計したディジタルフィルタにおいて,階段関数の信号によって連続時間信号を近似した事により生じる誤差の補正方法を示す.ここで,階段関数の信号で表示されるインパルス応答を,アナログフィルタのインパルス応答の近似と考える.この近似の精度は,アナログフィルタとディジタルフィルタにおける,時間 0 から ∞ の間のインパルス応答の面積比から与えられる.

アナログフィルタの伝達関数G (s),そのインパルス応答を g (t),ディジタル化されたフィルタのシステム伝達関数H (z),そのインパルス応答の離散時間信号を h [n] とする.ここで,時間 0 から ∞ におけるインパルス応答 g (t) の面積は以下で与えられる.

(10-5) 0gtdt=limt0tgτdτ=lims0Gs

上記は一見不思議な式であるが,表5-2 に示したラプラス変換積分公式と最終値定理から以下の通り求められる.

(10-6) L0tfτdτ=1sFs limtft=lims0sFs limt0tfτdτ=lims0s1sFs=lims0Fs

一方,10.2 に基づいてゲインが補正されたシステム伝達関数 T H (z) における,インパルス応答の離散時間信号は T h [n] となる.このゲインが補正されたインパルス応答のインパルス列の面積は,階段関数の信号によるインパルス応答の面積と等しいため,時間 0 から ∞ における階段関数の信号によるインパルス応答の面積は以下で与えられる.

(10-7) n=0Thn=Tlimz1Hz

従って,階段関数の信号による誤差を補正したシステム伝達関数 HN (z) は以下で与えられる.

(10-8) HNz=lims0GsTlimz1HzTHz=lims0Gslimz1HzHz

Note: 本資料では以降,実装モデルにおいてゲインが適正となる様に補正したシステム伝達関数を,正規化システム伝達関数 HN (z) と記す.

10.3 の伝達関数の補正において周波数伝達関数の直流の際のゲインを論じたが,これは伝達関数やボード線図を読み解く際の技法である.一方,上記の正規化システム伝達関数における s → 0, z → 1 と,直流のゲインを求める際の ω → 0 は同じ事を意味している.即ち,ディジタルフィルタのインパルス応答が,アナログフィルタのインパルス応答を近似する様にシステム伝達関数を補正する方法は,これらの直流におけるゲインが等しくなる様に補正する技法と論理的に等価となる.この方法は周波数伝達関数の直流におけるゲインが 0 倍となるか,発散し無い限り適用できる.

Note: 式 (10‑4) は,正規化された周波数伝達関数における直流のゲインを示しているが,以上の議論より,これは正規化システム伝達関数 HN (z) のインパルス応答を表示する階段関数の信号の面積と等しい.従って,階段関数の信号による連続時間信号の近似が,この誤差の主な要因となる.

10.5. 誤差を補正したローパスフィルタの例

9.10 に述べたローパスフィルタのシステム伝達関数を,上記のアナログフィルタとディジタルフィルタのインパルス応答の関係に基づいて補正した例を示す.9.10 に述べたローパスフィルタの伝達関数を以下に再掲する.

(10-9) Gs=1τs+1=ωcs+ωc

また,上記の伝達関数に対応するシステム伝達関数を以下に再掲する.

(10-10) Hz=ωc1-e-ωcTz-1

従って,正規化システム伝達関数は以下で与えられる.

(10-11) HNz=lims0Gslimz1HzHz=1-e-ωcT1-e-ωcTz-1

上記の正規化システム伝達関数に対応する,周波数伝達関数を以下に示す.

(10-12) HNejωT=1-e-ωcT1-e-ωcTe-jωT

上記の周波数伝達関数のゲイン線図を図10‑6に示す.横軸の角周波数は ω / ωc として正規化されている.また図中の青線は正規化されたナイキスト周波数 fn を示しており,この例ではカットオフ周波数 ωc の 20 倍の角周波数としている.

f:id:suzumushi0:20170611122157p:plain
図10-6: 正規化されたローパスフィルタのゲイン線図の例.

上記のゲイン線図と図6-8 に示した一次遅れ要素 (ローパスフィルタ) のゲイン線図を比較すると,ナイキスト周波数以下の特性は高い精度で一致する.

10.6. 0 次ホールドの伝達関数

ディジタル信号処理の論理モデルにおける 0 次ホールドは,インパルス列を階段関数の信号に変換する.単位インパルス関数は,単位ステップ関数の微分であるから,0 次ホールドはインパルス列中の単位インパルス関数を積分してステップ関数とし,サンプリング周期 T が経過した時点でそれを反転する事によってインパルス列を階段関数の信号に変換する.従って 0 次ホールドの伝達関数 G (s) は以下で与えられる.

(10-13) Gs=1s-e-sTs=1-e-sTs

10.2 に述べた通り,上記の正規化伝達関数 GN (s) は以下で与えられる.

(10-14) GNs=1TGs=1-e-sTsT

従って,上記の周波数伝達関数は以下の通りとなる.

(10-15) GNjω=1-e-jωTjωT

この 0 次ホールドの伝達関数は,システム伝達関数や正規化システム伝達関数には含まれていないため,ディジタルフィルタを設計する際には,正規化システム伝達関数に基づく周波数伝達関数に加えて,0 次ホールドの周波数伝達関数も考慮しなければならない事に注意せよ.

Note: 0 次ホールドの伝達関数を D/A 変換器の伝達関数と誤解しない事.D/A 変換器は離散時間信号を階段関数の信号に変換するが,これらの関係は定数係数線形常微分方程式で記述出来ないため,伝達関数も周波数伝達関数も定義されない.

上記の周波数伝達関数ボード線図を図10‑7に示す.横軸の角周波数は ωT として正規化されている.また図中の青線は正規化されたナイキスト周波数を示している.

f:id:suzumushi0:20170611122158p:plain
図10-7: 0 次ホールドのボード線図.

図から明らかな通り,0 次ホールドのゲインは,ナイキスト周波数において約 −3.9 [dB] (約 0.64倍),位相は −90 [deg] となる.この様なゲインの低下,位相の遅れを避けるために,D/A 変換器のサンプリング周波数を,ディジタル信号処理のサンプリング周波数の数倍として,更に離散時間信号を補間して D/A 変換器の入力とする技法が用いられる場合がある.

本章の参考文献

筆者, "修士論文," 某大学大学院工学研究科, 1986.

Note: 筆者の修論のテーマは自動制御理論でもディジタル信号処理でも無い.無論これらを道具として最大限活用したが,研究対象とはしていない.10.3, 10.4, 10.5は,実験データを分析する過程で,データをフィルタリングすると精度が低下する事象が生じたため,一晩徹夜して問題を解決した際に見出したものである.

9. ディジタル信号処理と Z 変換

本章では,先ず9.1において,サンプリング定理を復習する.次に,9.2において,ディジタル信号処理におけるサンプラーの論理モデルとサンプル値の表現を,9.3 において Z 変換の導出を示す.そして 9.4 以降において Z 変換の定義,公式と定理,ラプラス変換の像関数からの Z 変換,差分方程式の解法,逆 Z 変換,システム伝達関数と周波数応答を復習する.これらは教科書に十分な説明が記載されているため,本資料では数式の証明や導出等は省略し,留意すべき点を Note に示す.最後に 9.10 において,以上のディジタル信号処理の教科書的説明に基づいて設計したローパスフィルタの例を述べ,教科書通りの設計では期待した特性のフィルタは得られない事を示す.

9.1. サンプリング定理

サンプリング (標本化) とは連続時間信号を離散時間信号に変換する方法の事である.多くの場合,連続時間信号の瞬時値を周期的に計測して,これを離散時間信号とする.離散時間信号を計測する周期はサンプリング周期,その逆数はサンプリング周波数と呼ばれる.

変換された離散時間信号から,連続時間信号に含まれる f0 までの周波数成分を再現するには,サンプリング周波数 fsfs ≥ 2 f0 でなければならない.サンプリング周波数 fs の 1/2 の周波数をナイキスト周波数と言う.

離散時間信号では,連続時間信号に含まれるナイキスト周波数 fn = fs / 2 以下の周波数成分 f0 と,ナイキスト周波数を超える周波数成分 N fs ± f0, (N = 1, 2, 3, ...) を区別する事はできない.これは,異名現象 (エイリアシング) と呼ばれる.従って,連続時間信号をサンプリングする場合は,サンプラーの前にアナログフィルタ等を設けて,ナイキスト周波数を超える周波数成分を十分に減衰させる,若しくはナイキスト周波数を超える周波数成分によって生じるノイズを許容する事が必要となる.

9.2. サンプル値の表現

ディジタル信号処理やサンプル値制御におけるサンプラーの論理モデルを図9‑1 に示す.論理モデルにおけるサンプラーは,連続時間信号の入力と,サンプリング周期 T = 1 / fs で生成される単位インパルス関数 (櫛型関数) を乗算する事によって,連続時間信号の入力の瞬時値をインパルス列として抽出する.

f:id:suzumushi0:20170611122150p:plain
図9-1: サンプラーの論理モデル.

即ち,連続時間信号の入力を x (t), (t ≥ 0),インパルス列を p (t), (t ≥ 0) とすると,p (t) は以下の式で記述される.ここで,δ (t) は,5.2 に述べた単位インパルス関数 (デルタ関数) である.

(9-1) pt =xtδt+δt-T+δt-2T+δt-3T+ =xtn=0δt-nT =n=0xnTδt-nT,t0

Note: 単位インパルス関数は理論上,想像上の関数であるため,サンプラーの論理モデルと現実の実装の間には相当な乖離がある事に留意せよ.また,この論理モデルは,連続時間信号の入力 x (t) と周期 T で生成される単位インパルス関数を乗算しているため,連続時間信号の入力 x (t) とインパルス列 p (t) の関係は定数係数線形常微分方程式では記述できない.従って,これらの間には伝達関数や周波数伝達関数は定義されない事にも注意せよ.

9.3. Z 変換の導出

次に,連続時間信号の入力とインパルス列の関係から,Z 変換を導出する.式 (9-1) のラプラス変換を以下に示す.変換には,表5-2 に示した時間軸上の平行移動定理を利用している.

(9-2) Ps=n=0xnTe-nTs

上記の式において z = eTs に置き換えると,以下の通り Z 変換が得られる.

(9-3) Pz=n=0xnTz-n

ここで,連続時間信号の入力 x (t) に対して,x [n] = x (nT) となる数列を離散時間信号とする.x [n] は離散時間信号の n 番目の標本となる.この場合の p (t) 及び P (z) を以下に示す.

(9-4) pt=n=0xnδt-nT,t0 Pz=n=0xnz-n

以上から明らかな通り,Z 変換とはインパルス列のラプラス変換の別表現である.また,インパルス列における各単位インパルス関数の係数が,離散時間信号となる事に注意せよ.

9.4. Z 変換の定義

n ≥ 0 で定義された数列 (離散時間信号) f [n] の Z 変換 F (z) を以下に示す.

(9-5) Fz=n=0fnz-n

例えば f [0] = 3, f [1] = 1, f [2] = 4, f [n] = 0, (n > 2) となる数列 f [n] の Z 変換は以下の通りとなる.

(9-6) Fz=n=0fnz-n=3+z-1+4z-2

ここで z 複素数z = eTs = eT(σ + ) = re jωT と表す事が多い.F (z) が収束する z の範囲を収束域と言う.f [n] と F (z) の関係を以下の様に表示する場合もある.

(9-7) Zfn=Fz

同様に,f [n] は F (z) の逆 Z 逆変換とも呼ばれ,以下の様に表示される.

(9-8) Z-1Fz=fn

逆 Z 変換の公式を以下に示す.ここで,CF (z) の収束領域内にあり,原点の周りを反時計方向に回る閉路である.

(9-9) fn=1j2πCFzzn-1dz

Z 変換は差分方程式の解法やディジタル信号処理の解析に用いられるが,フーリエ解析の様に観測した離散時間信号を Z 変換する事は稀である.このため,ラプラス変換と同様に,上記の逆 Z 変換の公式を実際に解く必要は殆ど無く,多くの場合逆変換には以下の Z 変換の公式や定理が用いられる.

9.5. Z 変換の公式と定理

Z 変換の主要な公式を以下に示す.

表9-1: Z変換の主要な公式.
項番 fn,n0 Fz 収束域
1 δn=1,n=00,n0 1  
2 un=1,n00,n<0 11-z-1 z>1
3 1 11-z-1 z>1
4 n z-11-z-12 z>1
5 αn 11-αz-1 z>|α|
6 nαn αz-11-αz-12 z>|α|

項番 1 は単位インパルス数列と呼ばれる.あの難解な単位インパルス関数が離散時間信号ではこんなに簡単になると喜んでいる者は,9.2 からの説明を再度読む事.離散時間信号とはインパルス列における各単位インパルス関数の係数である.項番 1 の数列を係数とするインパルス列は,単位インパルス関数となるため,これを単位インパルス数列と称しているだけの事である.

項番 2 は単位ステップ数列と呼ばれる.f [n] の定義域は n ≥ 0 であるから,これを項番 3 の様に表す場合もある.また,項番 5, 6 の α には複素数が含まれる.尚,項番 2 以降は,何れも高校数学で学習した以下の等比数列の和の公式から求める事ができる.

(9-10) n=0arn=a1-r,r<1

最後に,Z 変換の主要な定理を以下に示す.

表9-2: Z変換の主要な定理.
線 形 性 Zafn+bgn=aZfn+bZgn=aFz+bGz
畳み込み Zk=0fn-kgk=ZfnZgn=FzGz
時間軸のシフト Zfn-k=z-kZfn=z-kFz,k0
指数数列の乗算 Zanfn=Fa-1z
F (z) の微分 Znfn=-zddzFz
終値定理 f=limz11-z-1Fz
初期値定理 f0=limzFz

9.6. ラプラス変換の像関数からの Z 変換

ラプラス変換の原関数とその像関数,及び原関数を周期 T でサンプリングした数列とその Z 変換の対応関係を表9‑3 に示す.この表から,時間関数を周期 T でサンプリングした数列の Z 変換を,その時間関数のラプラス変換から求める事ができる.

表9 3: ラプラス変換と Z 変換の対応関係.
原関数 ft , t 0 像関数 Fs fn,n0 Fz
ut=1,t00,t<0 1s un=1,n00,n<0 11-z-1
t 1s2 nT Tz-11-z-12
e-at 1s+a e-naT 11-e-aTz-1

例として,連続時間関数 eat cos (ωt) を周期 T でサンプリングした数列の Z 変換を,以下のラプラス変換の像関数から求める.

(9-11) Le-atcosωt=s+as+a2+ω2

ここで,右辺は以下の通り部分分数に分解できる.

(9-12) s+as+a2+ω2=121s+a+jω+1s+a-jω

表9-3 より,上記の像関数に対応する Z 変換を以下に示す.

(9-13) Zs+as+a2+ω2 =1211-e-a+jωTz-1+11-e-a-jωTz-1 =122-e-a+jωT+e-a-jωTz-11-e-a+jωTz-11-e-a-jωTz-1 =122-e-aTe-jωT+ejωTz-11-e-aTe-jωT+ejωTz-1+e-2aTz-2 =1-e-aTcosωTz-11-2e-aTcosωTz-1+e-2aTz-2

別の例として,連続時間関数 eat sin (ωt) のラプラス変換を以下に示す.

(9-14) Le-atsinωt=ωs+a2+ω2

上記の例と同様に,これを周期 T でサンプリングした数列の Z 変換は以下の通りとなる.

(9-15) Zωs+a2+ω2=e-aTsinωTz-11-2e-aTcosωTz-1+e-2aTz-2

Note: 原関数が単位インパルス関数となる場合は,これをサンプリングしても離散時間信号は無限大に発散してしまうため,Z 変換の単位インパルス数列には対応しない事に注意せよ.

9.7. Z 変換による差分方程式の解法

Z 変換による下記の 1 階の差分方程式 (定数係数線形差分方程式) の解法を示す.高階の場合も解法は原理的に同じである.

(9-16) yn-ayn-1=xn

上記の両辺の Z 変換を以下に示す.

(9-17) Yz-az-1Yz=Xz

ここで,X (z), Y (z) は各々数列 x [n], y [n] の Z 変換を表す.上記を Y (z) について解くと以下の通りとなる.

(9-18) Yz=Xz1-az-1

上記は数列 y [n] の Z 変換を与えるから,これを下記の様に逆 Z 変換する事によって y [n] が得られる.

(9-19) yn=Z-1Xz1-az-1

Note: 式 (9‑16) において,y [0] を求める際に y [−1] が必要となるが,ディジタル信号処理では,n < 0 における離散時間信号 f [n] は 0,即ち初期条件を 0 とする事を暗黙の前提としている.これは,表9‑2 の時間軸シフトの定理は f [n] = 0, (n < 0) を前提としている事からも明らかである.ディジタル信号処理ではまず必要とされないが,差分方程式の初期値問題を解く場合は,表9‑2 の時間軸シフトの定理を k < 0 とした場合が必要となる.参考までに k = −2 の場合の時間軸シフトの定理を以下に示す.

(9-20) Zfn+2=z2Zfn-z2f0-zf1=z2Fz-z2f0-zf1

9.8. 逆 Z 変換

逆 Z 変換の具体的な方法を説明する.説明を簡単にするために,Y (z) がべき級数の比で与えられる下記の関数の逆変換方法を示す.

(9-21) Yz=b0+b1z-1a0+a1z-1+a2z-2

べき級数展開による方法

Y (z) がべき級数の比で与えられる場合は,以下の様に,分子を分母で除算するべき級数展開によって,逆 Z 変換を求める事ができる.

(9-22) Yz =b0+b1z-1a0+a1z-1+a2z-2 =y0+y1z-1+y2z-2+

部分分数分解による方法

ラプラス変換と同様に,逆 Z 変換においても,式 (9‑21) を部分分数に分解して,Z 変換の公式に当てはめる事によって簡易に逆変換ができる.Y (z) の極が単極,若しくは複素数の場合,これらを各々 α, β とすると,式 (9‑21) は以下の様に部分分数に分解できる.

(9-23) Yz=k0+k11-αz-1+k21-βz-1

ここで,k1, k2 は以下の様に留数定理 (ヘビサイドの展開定理) から,k0 は係数比較から求められる.

(9-24) k1=limzα1-αz-1Yz=1-αz-1Yzz=α k2=limzβ1-βz-1Yz=1-βz-1Yzz=β k0=b0a0-k1+k2

以上と Z 変換表より,y [n] は以下の通りとなる.

(9-25) yn=k0δn+k1αn+k2βn

尚,Y (z) の極が n 重極の場合,これを α とすると,Y (z) は以下の様に部分分数に分解できる.

(9-26) Yz=k11-αz-1n+k21-αz-1n-1++kn1-αz-1

ここで,kr は以下の留数定理から求められる.

(9-27) kr =1r-1!limzαdr-1dzr-11-αz-1nYz =1r-1!dr-1dzr-11-αz-1nYzz=α

9.9. システム伝達関数と周波数応答

9.3 に述べた通り,Z 変換とはインパルス列のラプラス変換の別表現であり,かつ離散時間信号の初期条件は 0 であるから,システムに入出力されるインパルス列のラプラス変換の比,即ち Z 変換の比から伝達関数や周波数伝達関数が定義できる.システムの入力の離散時間信号 x [n] と出力の離散時間信号 y [n] が以下の差分方程式で記述されるシステムにおける,入力 x [n] と出力 y [n] の Z 変換の比を,システム伝達関数と言う.

(9-28) yn=k=1Nakyn-k+r=0Mbrxn-r

上記の差分方程式は,出力の離散時間信号 y [n] が,過去 N 回の出力と,現在から過去 M 回の入力の一次結合となるシステムを示している.上記の Z 変換を以下に示す.

(9-29) Yz=k=1NakYzz-k+r=0MbrXzz-r

ここで,X (z), Y(z) は各々入出力の離散時間信号 x [n], y [n] の Z 変換を表す.よってシステム伝達関数 H (z) は以下で与えられる.

(9-30) Hz=YzXz=r=0Mbrz-r1-k=1Nakz-k

また,Z 変換はラプラス変換において z = eTs に置き換えたものであるから,周波数伝達関数は,以下に示す様に,システム伝達関数において z = e jωT に置き換える事によって得られる.

(9-31) HejωT=r=0Mbre-jrωT1-k=1Nake-jkωT

9.10. ローパスフィルタの例

以上,教科書的なディジタル信号処理の基本を復習した.ここで,ディジタルフィルタの設計の各論に移る前に,以下の伝達関数で示される,時定数が τ,カットオフ周波数が ωc = 1 / τ のローパスフィルタをディジタル化した例を示す.

(9-32) Gs=1τs+1=ωcs+ωc

表9‑3 より,上記の伝達関数に対応するシステム伝達関数を以下に示す.サンプリング周期は T とする.

(9-33) Hz=ωc1-e-ωcTz-1

また,上記のシステム伝達関数に対応する,周波数伝達関数を以下に示す.

(9-34) HejωT=ωc1-e-ωcTe-jωT

上記の周波数伝達関数のゲイン線図を図9‑2 に示す.横軸の角周波数は ω / ωc として正規化されている.また図中の青線は正規化されたナイキスト周波数 fn を示しており,この例ではカットオフ周波数 ωc の 20 倍の角周波数としている.尚,サンプリング周波数は fs = 1 / T = 1 [MHz] としている.

f:id:suzumushi0:20170611122151p:plain
図9-2: ディジタル化されたローパスフィルタのゲイン線図の例.

上記のゲイン線図と図6-8 に示した一次遅れ要素 (ローパスフィルタ) のゲイン線図を比較すると,ナイキスト周波数以下の特性は一見似ている様に見えるが,ゲインの最大値が 120 dB 以上となっており特性が全く異なるローパスフィルタとなっている.

実は,ディジタル信号処理の論理モデルは,現実の実装モデルと相当な乖離があるため,理論モデルを解説した教科書通りに設計しても,期待した特性のディジタルフィルタは得られない.そこで,10章において,論理モデルと実装モデルの関係を述べ,論理モデルを実装モデルに適用する方法を明らかにする.

10. ディジタル信号処理の論理モデルと実装モデル - 趣味人のブログ

本章の参考文献

明石一, "制御工学 増訂版," 大学講座 機械工学12, 共立出版, 1979.
A. Oppenheim, R. Schafer, 伊達訳, "ディジタル信号処理," (上), (下), コロナ社, 1978.

8. フーリエ級数展開とフーリエ変換 (続き)

参考: フーリエ変換と単位インパルス関数

8.2 に示したフーリエ変換の定義では,ごく基本的な関数でもフーリエ変換が収束しない場合がある.例えば f (t) = 1 や周期関数の f (t) の変換が収束しない事は,定義式 (8‑19) から自明である.しかし,5.2に示した単位インパルス関数 (デルタ関数) によって,これらの関数のフーリエ変換を形式的に表示する事が出来る.

単位インパルス関数のフーリエ変換

式 (5‑5) に示した単位インパルス関数 δ (t) の定義より,これのフーリエ変換は以下で与えられる.

(8-26) Fδt=-δte-jωtdt=e0=1

上記は,時間関数の単位インパルス関数が,定数の周波数関数,即ち振幅が等しく周波数成分が ∞ までとなる正弦波と直流成分の重ね合わせによって表せる事を示している.上記のフーリエ逆変換は以下で与えらえる.

(8-27) δt=12π-ejωtdω

ここで,上記の式の変数を入れ替えると,以下の通り変形できる.

(8-28) δω=12π-e-jωtdt

式 (8-27), (8-28) の積分は収束しないが,単位インパルス関数はこれらの式によって形式的に表示する事ができる.これらの怪しさ爆発する式は,単位インパルス関数の積分表示と呼ばれる.

直流成分のフーリエ変換

単位インパルス関数の積分表示より f (t) = Aフーリエ変換は以下で表示する事ができる.

(8-29) Fω=-Ae-jωtdt=2πAδω

上記のフーリエ逆変換を以下に示す.

(8-30) F-1Fω=12π-2πAδωejωtdω=Ae0=A

上記は,定数の時間関数,即ち直流成分が,周波数関数の単位インパルス関数によって表せる事を示している.これは前に述べた,時間関数の単位インパルス関数が,定数の周波数関数によって表せる事と対称的な関係を示している.

正弦関数,余弦関数のフーリエ変換

単位インパルス関数の積分表示より,f (t) = A sin (ω0t) のフーリエ変換は以下で表示する事ができる.

(8-31) Fω =-Asinω0te-jωtdt =A-ejω0t-e-jω0tj2e-jωtdt =πAj-e-jω-ω0t-e-jω+ω0t2πdt =πAjδω-ω0-δω+ω0 =jπAδω+ω0-δω-ω0

上記のフーリエ逆変換を以下に示す.

(8-32) F-1Fω =12π-jπAδω+ω0-δω-ω0ejωtdω =jA2-δω+ω0ejωtdω--δω-ω0ejωtdω =Aejω0t-e-jω0tj2 =Asinω0t

同様に,単位インパルス関数の積分表示より,f (t) = A cos (ω0t) のフーリエ変換は以下で表示する事ができる.

(8-33) Fω =-Acosω0te-jωtdt =A-ejω0t+e-jω0t2e-jωtdt =πA-e-jω-ω0t+e-jω+ω0t2πdt =πAδω-ω0+δω+ω0

ここで,上記の F (ω0) は以下で与えられる.

(8-34) Fω0=πaδ0

上記及び,直流成分,正弦関数,余弦関数のフーリエ変換を表示した式 (8-29), (8-31), (8-33) からも明らかな通り,これらのフーリエ変換には単位インパルス関数が現れる.従って,8.2 に述べた方法では,これらのフーリエ変換から正弦波の振幅や初期位相,直流成分を求める事はできない.勿論,ご都合主義的に,これらの式の単位インパルス関数を無視する事も可能であるが,フーリエ変換の結果に単位インパルス関数を含む項と含まない項が混在した場合の解釈が困難となる.

また,これらの結果に対してパーセバルの定理は成立しないとされている.そもそも,理論上,想像上の関数である単位インパルス関数のエネルギーを論じる事の意味を見い出す事は困難である.結局,単位インパルス関数は理論上,想像上の関数であるため,これが時間関数として現れても,周波数関数として現れても,その意味を正確に議論する事は容易ではない.従って,単位インパルス関数によって,直流成分や周期関数のフーリエ変換が「形式的に表示できるだけ」と考えた方が現実的である.

では,何故この様な無用と思われる議論を述べているかと言うと,9.2 に述べるサンプル値の表現にはこの表示方法が使われている点と,以下の様な議論が可能となる点にある.

単位ステップ関数のフーリエ変換

以下で定義される単位ステップ関数 u (t) のフーリエ変換を求める.

(8-35) ut=1,t10,t<0

式 (5-8) に示した通り,単位インパルス関数は単位ステップ関数の微分であるから,表8-1 に示したフーリエ変換微分の定理により,単位ステップ関数のフーリエ変換は以下で表示する事ができる.

(8-36) jωFut=Fddtut=Fδt=1 Fut=1jω

以上が,教科書的な単位ステップ関数のフーリエ変換の説明であるが,残念な事に単位ステップ関数は不連続な関数であるから,上記をフーリエ逆変換すると不連続点の u (0) は近似となり,単位ステップ関数に厳密には収束しない.前に述べたフーリエ逆変換の収束の規則から,上記をフーリエ逆変換した関数 uc (t) を以下に示す.

(8-37) uct=1,t>012,t=00,t<0

次に,改めて上記の関数 uc (t) のフーリエ変換を求める.これには以下の符号関数 sgn (t) を用いる.

(8-38) sgnt=1,t>00,t=0-1,t<0

符号関数 sgn (t) のフーリエ変換は以下の通りとなる事が既知である.(以下の導出は工学の範囲を超えるため割愛する).

(8-39) Fsgnt=2jω

ここで,関数 uc (t) と符号関数 sgn (t) には以下の関係がある.

(8-40) uct=sgnt+12

よって,フーリエ変換の線形性の定理と,前に述べた f (t) = aフーリエ変換から,関数 uc (t) のフーリエ変換は以下で表示する事ができる.

(8-41) Fuct=Fsgnt2+F12=1jω+πδω

上記の πδ (ω) は 1 / 2 の直流成分を表す.従って,式 (8‑36) に示した,単位ステップ関数 u (t) のフーリエ変換には,単位ステップ関数の直流成分が正確に表示されていない事が判る.

余談: ラプラス変換 = フーリエ変換

一部の教科書に,ラプラス変換において s = とすればフーリエ変換となるとの説が見受けられる.おおむね,この説はラプラス変換積分範囲を −∞ から ∞ とした両側ラプラス変換に基づいている.ラプラス変換の原関数 f (t) の両側ラプラス変換 F (s) の定義を以下に示す.式 (5‑1) に定義したラプラス変換と下記を区別するために,ここでは前者を片側ラプラス変換と呼ぶ.

(8-42) Bft=Fs=-fte-stdt

この説では,上記の両側ラプラス変換において s = とすればフーリエ変換の定義式となり,かつ原関数を f (t) = 0, (t < 0) とすれば片側ラプラス変換の定義式と等価となるから,片側ラプラス変換において s = とすればフーリエ変換となるとしている.(尚,ここまで丁寧に書かれていない説明が殆どである).

片側ラプラス変換と両側ラプラス変換の定義は似ているため,表5‑1や表5‑2に示した片側ラプラス変換の公式や定理の殆どは両側ラプラス変換においても成立するが,残念ながら成立しない場合もある.定理が異なる例として片側ラプラス変換における微分の定理を以下に示す.

(8-43) Lddtft=sFs-f+0

一方,両側ラプラス変換における微分の定理は以下となる.

(8-44) Bddtft=sFs

当然,フーリエ変換微分の定理も上記と同様となる.

(8-45) Fddtft=jωFω

これらから明らかな通り,片側ラプラス変換と異なり,両側ラプラス変換フーリエ変換における微分の定理には − f (+0) の項が無いため,原関数を f (t) = 0, (t < 0) とする制約は,f (t) の片側ラプラス変換と両側ラプラス変換が等価となる条件としては不十分であり,これに加えて少なくとも f (0) = 0 とする制約が必要となる.即ち,両側ラプラス変換は,片側ラプラス変換と似てはいるが,これを包含するものでは無い

5.3に述べた通り,片側ラプラス変換によって定数係数線形常微分方程式の初期値問題における特殊解を求める事ができる.この場合,原関数 f (t) の定義域は t ≥ 0 となる制約があるが,任意の初期条件を与える事ができる.一方,4.4に述べた未定係数法と同様に,フーリエ変換によって定数係数線形常微分方程式の非同次方程式の特殊解 (定常解)  yp (t) を求める事ができる (これは実質的に未定係数法と同じ解法であり,わざわざフーリエ変換によって解く意味は無いため詳細は割愛する).この場合,未定係数法と同様に,原関数 f (t) の定義域に制約は無く解は初期条件に依存しない.

従って,片側ラプラス変換フーリエ変換を定数係数線形常微分方程式の解法とした場合,これらは異なる制約条件における異なる問題の解法となるため,単に原関数を f (t) = 0, (t < 0) とし s = とすればこれらは等価という説明は著しく正確性を欠く

多分この説の出処は,伝達関数において s = とすれば周波数伝達関数となる事を,ラプラス変換フーリエ変換の関係として説明したい点にあると推測するが,伝達関数は初期条件が全て 0,即ち f (0) = f ' (0) = ... = 0 を前提としており,この場合は片側ラプラス変換と両側ラプラス変換,即ちフーリエ変換は等価となるため,伝達関数において s = とすれば周波数伝達関数となるのである

8. フーリエ級数展開とフーリエ変換

本章はフーリエ級数展開フーリエ変換について述べる.これらフーリエ解析は振動の分析に適用する事ができ,線形時不変システムより広い分野で利用されているため,非常に多くの教科書がある.このため,本資料では数式の証明や導出は省略し,振動の分析という観点からの解説を行う.尚,フーリエ解析における種々の関係は,以下の f (t) が複素数の場合でも成立するが,その物理的な意味を見い出す事は容易ではないため,本資料では扱わない.

8.1. フーリエ級数展開

周期 2π の周期関数のフーリエ級数展開

f (t) を区間 [−π, π] で定義された周期 2π の周期関数とすると,以下の式に示す通り,f (t) は定数,及び周期 2π / n, (n = 1, 2, 3, ...) の正弦関数,及び余弦関数の和から成るフーリエ級数表す事ができる.下記の式の ~ は,フーリエ級数が厳密に f (t) に収束する場合と,近似となる場合がある事を示している.

(8-1) ft~a02+n=1ancosnt+bnsinnt

ここで,an 及び bnフーリエ係数と呼ばれ,以下で与えらえる.

(8-2) an=1π-ππftcosntdt,n=0,1,2, bn=1π-ππftsinntdt,n=1,2,3,

以上がフーリエ級数展開の教科書的定義である.ここで,式 (1-14) に示した三角関数の合成定理によって,式 (8-1) は以下の通り変形できる.

(8-3) ft~a02+n=1ancosnt+bnsinnt=a02+n=1Ansinnt+φn An=an2+bn2,φn=tan-1anbn,n=1,2,3,

即ち,f (t) を時間 [−π, π] で定義された周期 2π の振動における変位とすると,上記は角周波数が n,振幅が An,初期位相が φn, (n = 1, 2, 3, ...) となる正弦波と,直流成分 a0 / 2の重ね合わせで f (t) を表せる事を示している

式 (8‑2) のフーリエ係数を,以下の式に基づいて,複素フーリエ係数に変換すると,複素フーリエ級数が得られる.

(8-4) c0=a02,cn=an-jbn2,c-n= cn =an+jbn2

複素フーリエ級数と複素フーリエ係数を以下に示す.これらは,式 (8-1) 及び (8-2) の別表現となる.

(8-5) ft~n=-cnejnt cn=12π-ππfte-jntdt,n=0,±1,±2,

以上が,複素フーリエ級数展開の教科書的定義である.ここで,f (t) が実数の場合は cncn は共役複素数の関係となるため,式 (8‑5) は以下の通り変形できる.

(8-6) ft~n=-cnejnt =c0+n=1cnejnt+cne-jnt =c0+n=1j2cnejnt--j2cne-jntj2

ここで,An = | j2cn | = 2 | cn |, φn = Arg ( j2cn ) とすると,上記は更に以下の通り変形できる.

(8-7) ft~n=-cnejnt =c0+n=1Anejφnejnt-e-jφne-jntj2 =c0+n=1Anejnt+φn-e-jnt+φnj2 An=j2cn=2cn,φn=Argj2cn,n=1,2,3, cn=12π-ππfte-jntdt,n=0,1,2,

上記の複素フーリエ級数の総和における各項は,式 (2‑27) に示した複素指数関数による正弦波の記述となる.よって,f (t) を時間 [−π, π] で定義された周期 2π の振動における変位とすると,上記は角周波数が n,振幅が An,初期位相が φn, (n = 1, 2, 3, ...) となる正弦波と,直流成分 c0 の重ね合わせで f (t) を表せる事を示している.

Note: 上記の直流成分は c0 となるが,正弦波の振幅は 2 | cn | となる事に注意せよ.フーリエ級数に限らず,フーリエ解析を正弦波と直流成分の重ね合わせとして解釈する場合,直流成分と正弦波は区別して扱わなければならない

周期 2L の周期関数のフーリエ級数展開

f (t) を区間 [−L, L] で定義された周期 2L の周期関数とすると,以下の式に示す通り,f (t) は定数,及び周期 2L / n, (n = 1, 2, 3, ...) の正弦関数,及び余弦関数の和から成るフーリエ級数で表す事ができる.

(8-8) ft~a02+n=1ancosnπtL+bnsinnπtL

フーリエ係数 an 及び bn は以下で与えらえる.

(8-9) an=1L-LLftcosnπtLdt,n=0,1,2, bn=1L-LLftsinnπtLdt,n=1,2,3,

また,上記の別表現となる複素フーリエ級数と複素フーリエ係数を以下に示す.

(8-10) ft~n=-cnejnπtL cn=12L-LLfte-jnπtLdt,n=0,±1,±2,

以上の教科書的定義において L = π / ω と置き換えると,f (t) は区間 [−π / ω, π / ω] で定義された周期 2π / ω の周期関数となる.この場合,f (t) を時間 [−π / ω, π / ω] で定義された周期 2π / ω の振動における変位とすると,f (t) は,以下の通り,角周波数が ,振幅が An,初期位相が φn, (n = 1, 2, 3, ...) となる正弦波と,直流成分 a0 / 2の重ね合わせで表す事ができる.

(8-11) an=ωπ-π/ωπ/ωftcosnωtdt,n=0,1,2, bn=ωπ-π/ωπ/ωftsinnωtdt,n=1,2,3, An=an2+bn2,φn=tan-1anbn,n=1,2,3, ft~a02+n=1ancosnωt+bnsinnωt=a02+n=1Ansinnωt+φn

また,式 (8‑10) の f (t) が実数の場合は cncn は共役複素数の関係となるため,式 (8‑10) の複素フーリエ級数と複素フーリエ係数は以下の通り変形できる.

(8-12) cn=ω2π-π/ωπ/ωfte-jnωtdt An=2cn,φn=Argj2cn,n=1,2,3, ft~n=-cnejnωt=c0+n=1Anejnωt+φn-e-jnωt+φnj2

上記は角周波数が ,振幅が An,初期位相が φn, (n = 1, 2, 3, ...) となる正弦波と,直流成分 c0 の重ね合わせで f (t) を表せる事を示している.

フーリエ級数の収束

f (t) を区間 [a, b] で定義された周期関数とすると,f (t) が有限個の点を除いて連続であり,かつ全ての不連続点において,左側極限値と右側極限値が存在し,かつ a 及び b において各々左側極限値と右側極限値が存在する場合,f (t) を区分的に連続な関数と言う.また,周期関数 f (t) の導関数 f '(t) が区間 [a, b] において区分的に連続な関数となる場合,f (t) を区分的に滑らかな関数と言う.

区分的に連続,かつ区分的に滑らかな周期関数 f (t) のフーリエ級数の収束を以下に示す.

  1. 連続な点における f (t) のフーリエ級数f (t) に収束する.
  2. 不連続な点における f (t) のフーリエ級数は,下記の通り,不連続点の左側極限値と右側極限値の平均に収束する

(8-13) ft-0+ft+02

即ち,フーリエ級数は連続関数となる正弦波,及び直流成分の重ね合わせによって周期関数を表現するものであるため,連続な点は厳密に表現できるが,不連続点は近似となり,これを厳密に表現する事は出来ない

ギブスの現象

f (t) が以下で与えられる周期 2π の方形波 (矩形波) のフーリエ級数を例に,不連続点におけるフーリエ級数の収束過程を示す.

(8-14) ft=1,0t<π0,-πt<0

上記の方形波のフーリエ級数は以下で与えられる.

(8-15) ft~4πn=1sin2n-1t2n-1

上記のフーリエ級数の総和を n = 1 から m までの部分和で近似した際のグラフを以下に示す.

f:id:suzumushi0:20170611122147p:plain
図8-1: 方形波のフーリエ級数の収束.

図から明らかな通り,m が増大するに伴い,正弦波の重ね合わせは方形波に近づき,f (t) に収束してゆく.但し,不連続点の近傍で振動的な現象が観察される.m を更に増大させ,不連続点の近傍を拡大した図を以下に示す.

f:id:suzumushi0:20170611122148p:plain
図8-2: 方形波のフーリエ級数の収束(拡大).

図から明らかな通り,不連続点の近傍における振動的な現象は m が増大しても生じる.これはギブスの現象と呼ばれ,不連続点を含む周期関数のフーリエ級数の総和を部分和で近似した場合,不連続点の近傍で必ず生じる事が知られている.尚,不連続点 f (t) 近傍における,振動的な現象によるオーバーシュート,アンダーシュートの振幅は,不連続点の幅 f (t − 0) − f (t + 0) の約 0.09 倍となる.

以上が,ギブスの現象の教科書的説明である.ここで,振動の分析の観点からは,方形波のフーリエ級数の総和を部分和で近似する事は,方形波をローパスフィルタに通す事に等しい.従って,ギブスの現象は,カットオフ周波数を超える周波数成分を全く含まないローパスフィルタのインディシャル応答 (単位ステップ応答) や方形波の応答が,不連続点近傍で振動的となる事を示している.この振動はリップル (リンギング) と呼ばれる.

図6‑10 に示した通り,2 次遅れ要素のゲインは固有振動数 (共振周波数) 近傍で 0 [dB] 以上となるため,7.4 に示した通り,これのインディシャル応答が振動的となる事は直感的に理解できる.しかし,固有振動数を持たないローパスフィルタであっても,インディシャル応答が振動的となる場合がある事に十分留意する必要がある.事実,ローパスフィルタとして基本的な特性のバタワースフィルタは,次数を高くし減衰傾斜を急峻にすると,インディシャル応答や方形波の応答が振動的となり,リップルが生じる事が知られている.

リップルが生じるフィルタは,用途によっては許容されない場合がある.カットオフ周波数を超える周波数成分を全く含まないローパスフィルタは「理想的」フィルタと呼ばれ,またフィルタの減衰傾斜は急峻であるほど良い特性であるとの誤解が多く見受けられるが,過ぎたるは猶及ばざるが如し,その様な特性のフィルタは用途によっては問題を引き起こす可能性がある事に留意する必要がある.

窓関数によるリップルの抑制

リップルは,単純に減衰傾斜を緩やかにすれば抑制する事ができる.上に述べた方形波のフーリエ級数の総和を部分和で近似する場合は,部分和の各項に重み付の係数を乗じる事によってこれを実現できる.この様な重み付の技法は無数に存在するが,その代表的なものとして窓関数がある.最も良く使われる窓関数としてハミング窓を以下に示す.

(8-16) wn=0.54-0.46cos2πnN-1,0nN-1

ここで,w [n] は項 n の重み付の係数,N は項数を示す.上記は,多くの教科書に示されているハミング窓の定義で,n = N / 2 を中心とした項の重み付を高くする場合の定義式となる.

上に述べた,方形波のフーリエ級数の総和を部分和で近似する場合は,n = 0 を周辺の項の重み付を高くする必要があり,この場合のハミング窓の定義式は以下で与えられる.

(8-17) wn=0.54+0.46cosπnN-1,0nN-1

よって,ハミング窓で重み付された,方形波のフーリエ級数における m 項までの部分和による近似は以下で与えらえる.

(8-18) ft4πn=1m0.54+0.46cosπnm-1sin2n-1t2n-1

上記の計算結果のグラフを以下に示す.

f:id:suzumushi0:20170611122149p:plain
図8-3: ハミング窓で重み付された方形波のフーリエ級数の収束.

図から明らかな通り,ハミング窓による重み付により減衰傾斜が緩やかとなり,リップルは殆ど生じていない.その反面,重み付の影響により信号の立ち上がり,立下りに要する時間が長くなり,波形の変化が緩やかになっている.

8.2. フーリエ変換

フーリエ変換の定義

周期 2L の周期関数のフーリエ係数,フーリエ級数展開L→∞ とする事によってフーリエ変換フーリエ逆変換が得られる.区間 [−∞, ∞] で定義された関数 f (t) のフーリエ変換 F (ω) を以下に示す.

(8-19) Fω=-fte-jωtdt

関数 f (t) とフーリエ変換 F (ω) の関係を以下の様に表示する場合もある.

(8-20) Fft=Fω

同様に,関数 f (t) は F (ω) のフーリエ逆変換とも呼ばれ,以下の式で表す事ができるフーリエ級数と同様に,下記の式の ~ は,フーリエ逆変換が厳密に f (t) に収束する場合と,近似となる場合がある事を示している.

(8-21) ft~12π-Fωejωtdω

上記の関係を以下の様に表示する場合もある.

(8-22) F-1Fω~ft

フーリエ変換,逆変換の定義式は,係数 1 / (2π) を除けば,tω を入れ替えただけで,実質的に同じ式となる.

Note: フーリエ変換,逆変換の定義式は分野や文献によって異なる.特に,係数 1 / (2π) を変換,逆変換のどちらの係数とするかに差異が見られる.フーリエ係数との整合性の観点から,係数 1 / (2π) をフーリエ変換の定義式の係数とする場合もある.また,フーリエ変換,逆変換の定義式の対称性,パーセバルの定理によるエネルギーの計算との整合性の観点から,係数 1 / (2π) の平方根を変換,逆変換双方の係数とする場合もある.本資料ではラプラス変換の定義式との整合性の観点から,上記の様に係数 1 / (2π) を逆変換の定義式の係数としている.

以上が,フーリエ変換の教科書的定義である.ここで,f (t) が実数の場合は F (ω) と F (−ω) は共役複素数の関係となるため,フーリエ逆変換の定義式 (8‑21) は以下の通り変形できる.

(8-23) ft~12π-Fωejωtdω ft -F02π~12π0+Fωejωt+Fωe-jωtdω =1π0+jFωejωt--jFωe-jωtj2dω =0+Aωejωt+φω-e-jωt+φωj2dω Aω=Fωπ,φω=ArgjFωπ,ω>0

即ち,f (t) を時間 [−∞, ∞] で定義された振動における変位とすると,上記は角周波数をω,振幅を A (ω),初期位相を φ (ω), (ω > 0) とする正弦波と,直流成分 F (0) / (2π) の重ね合わせで f (t) を表せる事を示している.

Note: 以上述べた通り,フーリエ級数及びフーリエ変換は,ある波形の振動を正弦波の重ね合わせによって表現する際の,各正弦波の振幅,及び初期位相を示している.振動の分析にフーリエ解析を利用する際の目的によっては,角周波数と振幅の関係にしか興味が無い場合もあるが,特定分野の技術者向けの解説書ならいざ知らず,教科書レベルの文献において,初期位相は重要では無いと言った記述が散見されて仰天している.正弦波を合成した際に得られる波形は,各正弦波の振幅と初期位相に依存する.従って,振動を正弦波の重ね合わせによって表現する際に必要となる情報は,各正弦波の角周波数,振幅,及び初期位相となる.

フーリエ逆変換の収束

区間 [−∞, ∞] で定義された関数 f (t) が以下の式を満たす場合,f (t) を絶対積分可能な関数と言う.

(8-24) -ftdt<

関数 f (t) が絶対積分が可能,かつ区分的に連続,区分的に滑らかな場合,関数 f (t) のフーリエ変換 F (ω) が存在する.また,その場合のフーリエ逆変換の収束を以下に示す.

  1. 連続な点における f (t) のフーリエ逆変換は f (t) に収束する.
  2. 不連続な点における f (t) のフーリエ逆変換は,下記の通り,不連続点の左側極限値と右側極限値の平均に収束する

(8-25) F-1Fω=ft-0+ft+02

即ち,フーリエ級数と同様に,フーリエ逆変換も,連続な点は厳密に表現できるが,不連続点は元の関数の近似となり,これを厳密に表現する事は出来ない

フーリエ変換の定理

フーリエ変換の主要な定理を以下に示す.

表8-1: フーリエ変換の主要な定理.
線 形 性 Faft+bgt=aFft+bF[gt]=aFω+bGω
縮  尺 Ffat=1aFωa
対 称 性 FFt=2πf-ω
微  分 Fddtft=jωFω
畳み込み F-ft-τgτdτ=FωGω
パーセバルの定理 -ft2dt=12π-Fω2dω
時間軸のシフト Fft-a=e-jaωFft=e-jaωFω
周波数軸のシフト Fe-jatft=Fω+a

以下に続く.

8. フーリエ級数展開とフーリエ変換 (続き) - 趣味人のブログ