番外編: pdp11/70 レプリカ PiDP-11 (ソフトウエア)

前のブログを書いた直後に PiDP-11 のソフトが β 版となりましたので，これに含まれている UNIX について書きます．尚，全てソースコードが含まれて居ます．

UNIX 5, UNIX 6, UNIX 7

全く知らない世界です．まだ，コンソールがテレタイプ時代の OS の様です．

System III

大学，大学院時代はこれを使ってましたので懐かしいです．但し，私が使っていたマシンには vi や more が載ってましたが，これには含まれて居ません．多分，オリジナルリリースと商品版の差分なのでしょう．

System V

就職して最初に使った UNIX で，後に SystemV R1 と呼ばれるバージョンです．これにも pg (more の SysV 版) が含まれて居ません．尚，vi を動かす場合は ~/.profile に

TERM=vt100; export TERM

と書くと良いでしょう．

2.11 BSD

言わずと知れた Berkeley 版．TCP/IP も移植されており，インターネットアクセスできます．

番外編: pdp11/70 レプリカ PiDP-11 (ハードウエア)

PiDP-11 という，DEC 社製ミニコン pdp11/70 のレプリカの組み立てキットを購入しました．

pdp11/70 は Raspberry Pi 上のエミュレーター SimH で実現されます．PiDP-11 は，このエミュレータに接続されるコンソールとなります．

キットに含まれる部品．Raspberry Pi は含まれません．

基板に R，LED，ダイオード等を半田付けした状態．

基板にスイッチを半田付けした状態．スイッチの先頭が横一直線と成る様に，半田付けで調整する工程が結構大変．

基板裏に装着した Raspberry Pi 3 Model B+．

2.11 BSD が起動しました．

最初は，コンソールの SW からメモリに命令を書き込んで遊ぼうと思って居たのですが，RSX-11M や UNIX 等が動作するとあっては，そちらにのめり込んでしまいます．これは時間が幾らあっても足りませんｗ

13. FIR フィルタの概要

本章では，13.1 において FIR (Finite Impulse Response) フィルタ (有限インパルス応答フィルタ，非巡回型フィルタ) の基本を，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 変換を示す．

$$\begin{array}{}\text{(13-1)}& H\left(z\right)=\frac{Y\left(z\right)}{X\left(z\right)}=\sum _{n=0}^{N-1}h\left[n\right]{\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-n}\end{array}$$

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

$$\begin{array}{}\text{(13-2)}& y\left[n\right]=\sum _{k=0}^{N-1}h\left[k\right]\mathrm{<mspace\; width="0.1em"></mspace>}x\left[n-k\right]\end{array}$$

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

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

図13-1: FIR フィルタの基本構造．

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

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

$$\begin{array}{}\text{(13-3)}& H_d\left(e^{j\omega T}\right)=\sum _{n=0}^{\mathrm{\infty }}h\left[n\right]e^{-jn\omega T}\end{array}$$

ディジタル信号処理では，入出の離散時間信号 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) は以下で与えられる．

$$\begin{array}{}\text{(13-4)}& H\left(e^{j\omega T}\right)=\sum _{n=0}^{N-1}h\left[n\right]e^{-jn\omega T}\end{array}$$

ここで，インパルス応答を打ち切る事は，離散時間フーリエ変換 (DTFT) の総和 H_d (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] を以下に再掲する．

$$\begin{array}{}\text{(13-5)}& h\left[n\right]=\frac{1}{\tau }{e}^{-\frac{nT}{\tau }},\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{<mspace\; width="0.1em"></mspace>}\left(n\ge 0\right)\end{array}$$

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

$$\begin{array}{}\text{(13-6)}& H\left(z\right)=\sum _{n=0}^{N-1}h\left[n\right]z^{-n}=\frac{1}{\tau }\sum _{n=0}^{N-1}{e}^{-\frac{nT}{\tau }}{z}^{-n}={\omega }_c\sum _{n=0}^{N-1}{e}^{-n{\omega }_{c}T}{z}^{-n}\end{array}$$

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

$$\begin{array}{}\text{(13-7)}& H_N\left(z\right)=\frac{H\left(z\right)}{\underset{z\to 1}{\lim }H\left(z\right)}=\frac{1-e^{-{\omega }_{c}T}}{1-e^{-N{\omega }_{c}T}}\sum _{n=0}^{N-1}{e}^{-n{\omega }_{c}T}{z}^{-n}\end{array}$$

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

$$\begin{array}{}\text{(13-8)}& \begin{array}{}\begin{array}{r}{H}_N\left(e^{j\omega T}\right)\\ \end{array}\begin{array}{l}=\frac{1-e^{-{\omega }_{c}T}}{1-e^{-N{\omega }_{c}T}}\sum _{n=0}^{N-1}{e}^{-n{\omega }_{c}T}{e}^{-jn\omega T}\\ =\frac{\left(1-e^{-{\omega }_{c}T}\right)}{\left(1-e^{-N{\omega }_{c}T}\right)}\frac{\left\{1-e^{-N\left({\omega }_{c}T+j\omega T\right)}\right\}}{\left\{1-e^{-\left({\omega }_{c}T+j\omega T\right)}\right\}}\end{array}\end{array}\end{array}$$

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

図13-2: FIR によるローパスフィルタのボード線図の例．

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

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

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

$$\begin{array}{}\text{(13-9)}& w\left[n\right]=0.54+0.46\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\pi n}{N-1}\right),\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{<mspace\; width="0.1em"></mspace>}\left(0\le n\le N-1\right)\end{array}$$

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

$$\begin{array}{}\text{(13-10)}& h\left[n\right]=\frac{1}{\tau }\left\{0.54+0.46\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\pi n}{N-1}\right)\right\}{e}^{-\frac{nT}{\tau }},\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{<mspace\; width="0.1em"></mspace>}\left(n\ge 0\right)\end{array}$$

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

$$\begin{array}{}\text{(13-11)}& \begin{array}{}\begin{array}{r}H\left(z\right)\\ \end{array}\begin{array}{l}=\sum _{n=0}^{N-1}{h\left[n\right]z}^{-n}\\ ={\omega }_c\left[0.54\left(\sum _{n=0}^{N-1}{e}^{-n{\omega }_{c}T}{z}^{-n}\right)+0.23\sum _{n=0}^{N-1}\left\{e^{-n\left({\omega }_{c}T-\frac{j\pi }{N-1}\right)}+e^{-n\left({\omega }_{c}T+\frac{j\pi }{N-1}\right)}\right\}{z}^{-n}\right]\end{array}\end{array}\end{array}$$

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

$$\begin{array}{}\text{(13-12)}& \begin{array}{}\begin{array}{c}{H}_N\left(z\right)=\frac{H\left(z\right)}{\underset{z\to 1}{\lim }H\left(z\right)}\\ \underset{z\to 1}{\lim }H\left(z\right)={\omega }_c\left[0.54\frac{1-e^{-N{\omega }_{c}T}}{1-e^{-{\omega }_{c}T}}+0.23\left\{\frac{{1-e}^{-N\left({\omega }_{c}T-\frac{j\pi }{N-1}\right)}}{{1-e}^{-\left({\omega }_{c}T-\frac{j\pi }{N-1}\right)}}+\frac{{1-e}^{-N\left({\omega }_{c}T+\frac{j\pi }{N-1}\right)}}{{1-e}^{-\left({\omega }_{c}T+\frac{j\pi }{N-1}\right)}}\right\}\right]\end{array}\end{array}\end{array}$$

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

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

図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] の周波数成分の間には，周波数伝達関数の位相特性に応じた位相の遅れ進みが生じる．

$$\begin{array}{}\text{(13-14)}& y\left[n\right]=\sum _{k=0}^{\infty }h\left[k\right]\mathrm{<mspace\; width="0.1em"></mspace>}x\left[n-k\right]\end{array}$$

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

$$\begin{array}{}\text{(13-15)}& ry\left[n\right]=\sum _{k=0}^{\infty }h\left[k\right]\mathrm{<mspace\; width="0.1em"></mspace>}rx\left[n-k\right]=\sum _{k=0}^{\infty }h\left[k\right]\mathrm{<mspace\; width="0.1em"></mspace>}x\left[-n+k\right]\end{array}$$

次に，これらの離散時間信号 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] の関係は，以下の畳み込みで記述される．

$$\begin{array}{}\text{(13-16)}& rry\left[n\right]=ry\left[-n\right]=\sum _{k=0}^{\infty }h\left[k\right]\mathrm{<mspace\; width="0.1em"></mspace>}x\left[n+k\right]=\sum _{k=-\infty }^{0}rh\left[k\right]\mathrm{<mspace\; width="0.1em"></mspace>}x\left[n-k\right]\end{array}$$

上記の式から明らかな通り，散時間信号 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]，理論上，想像上のインパルス応答を h_d [n] = h_d [−n] とするフィルタの入出力の関係は以下の畳み込みで記述され，この入出力の離散時間信号間には位相の進みや遅れは生じず，即ち位相特性は角周波数によらず 0 となる．このため，想像上，理論上のインパルス応答が h_d [n] = h_d [−n] となるフィルタの位相特性は零位相特性となる．

$$\begin{array}{}\text{(13-17)}& y\left[n\right]=\sum _{k=-\infty }^{\infty }{h}_d\left[k\right]\mathrm{<mspace\; width="0.1em"></mspace>}x\left[n-k\right]\end{array}$$

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

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

$$\begin{array}{}\text{(13-18)}& h\left[n\right]=\left\{\begin{array}{l}{h}_d\left[n-\frac{N-1}{2}\right],\mathrm{<mspace\; width="0.2em"></mspace>}\left(0\le n\le N-1\right)\\ 0,\mathrm{<mspace\; width="6.8em"></mspace>}\left(n<0,\mathrm{<mspace\; width="0.2em"></mspace>}n\ge N\right)\end{array}\right.\end{array}$$

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

$$\begin{array}{}\text{(13-19)}& h\left[n\right]=h\left[N-1-n\right],\mathrm{<mspace\; width="0.2em"></mspace>}\left(0\le n\le N-1\right)\end{array}$$

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

図13-4: 線形位相 FIR フィルタの構造．

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

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

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

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

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

$$\begin{array}{}\text{(13-20)}& H_d\left(e^{j\omega T}\right)=\left\{\begin{array}{l}1,\mathrm{<mspace\; width="0.2em"></mspace>}\left(\left|\omega \right|\le {\omega }_c\right)\\ 0,\mathrm{<mspace\; width="0.2em"></mspace>}\left({\omega }_c<\left|\omega \right|\le \pi /T\right)\end{array}\right.\end{array}$$

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

$$\begin{array}{}\text{(13-21)}& \begin{array}{c}\begin{array}{r}{h}_d\left[n\right]\\ \\ \\ \end{array}\begin{array}{l}=\frac{T}{2\pi }{\int }_{-\frac{\pi }{T}}^{\frac{\pi }{T}}{H}_d\left(e^{j\omega T}\right)\mathrm{<mspace\; width="0.1em"></mspace>}{e}^{jn\omega T}d\omega \\ =\frac{T}{2\pi }{\int }_{-{\omega }_c}^{{\omega }_c}\mathrm{<mspace\; width="0.1em"></mspace>}{e}^{jn\omega T}d\omega \\ =\frac{1}{\pi n}\frac{e^{jn{\omega }_{c}T}-e^{-jn{\omega }_{c}T}}{j2}\\ =\frac{1}{\pi n}\mathrm{s}\mathrm{i}\mathrm{n}\left(n{\omega }_{c}T\right),\mathrm{<mspace\; width="0.2em"></mspace>}\left(n\ne 0\right)\end{array}\\ h_d\left[n\right]=\frac{T}{2\pi }{\int }_{-{\omega }_c}^{{\omega }_c}\mathrm{<mspace\; width="0.1em"></mspace>}{e}^{jn\omega T}d\omega =\frac{{\omega }_{c}T}{\pi },\mathrm{<mspace\; width="0.2em"></mspace>}\left(n=0\right)\end{array}\end{array}$$

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

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

$$\begin{array}{}\text{(13-22)}& h\left[n\right]=\left\{\begin{array}{l}\frac{\mathrm{s}\mathrm{i}\mathrm{n}\left\{\left(n-\frac{N-1}{2}\right){\omega }_{c}T\right\}}{\pi \left(n-\frac{N-1}{2}\right)},\mathrm{<mspace\; width="0.2em"></mspace>}\left(0\le n\le N-1,\mathrm{<mspace\; width="0.2em"></mspace>}n\ne \frac{N-1}{2}\right)\\ \frac{{\omega }_{c}T}{\pi },\mathrm{<mspace\; width="8em"></mspace>}\left(\mathrm{<mspace\; width="0.1em"></mspace>}n=\frac{N-1}{2}\right)\\ 0,\mathrm{<mspace\; width="10em"></mspace>}\left(n<0,\mathrm{<mspace\; width="0.2em"></mspace>}n\ge N\right)\end{array}\right.\end{array}$$

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

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

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

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

$$\begin{array}{}\text{(13-23)}& H\left(z\right)=\sum _{n=0}^{N-1}h\left[n\right]z^{-n}\end{array}$$

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

$$\begin{array}{}\text{(13-24)}& H\left(e^{j\omega T}\right)=\sum _{n=0}^{N-1}h\left[n\right]e^{-jn\omega T}\end{array}$$

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

図13-5: 線形位相 FIRフィルタのゲイン線図の例．

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

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

図13-6: 線形位相 FIRフィルタのインディシャル応答の例．

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

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

$$\begin{array}{}\text{(13-25)}& w\left[n\right]=0.54-0.46\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{2\pi n}{N-1}\right),\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{<mspace\; width="0.1em"></mspace>}\left(0\le n\le N-1\right)\end{array}$$

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

$$\begin{array}{}\text{(13-26)}& H\left(z\right)=\sum _{n=0}^{N-1}w\left[n\right]h\left[n\right]z^{-n}\end{array}$$

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

$$\begin{array}{}\text{(13-27)}& H_N\left(z\right)=\frac{H\left(z\right)}{\underset{z\to 1}{\lim }H\left(z\right)}=\frac{1}{\underset{z\to 1}{\lim }H\left(z\right)}\sum _{n=0}^{N-1}w\left[n\right]h\left[n\right]z^{-n}\end{array}$$

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

$$\begin{array}{}\text{(13-28)}& H_N\left(e^{j\omega T}\right)=\frac{1}{\underset{z\to 1}{\lim }H\left(z\right)}\sum _{n=0}^{N-1}w\left[n\right]h\left[n\right]e^{-jn\omega T}\end{array}$$

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

図13-7: 窓関数による線形位相 FIRフィルタのゲイン線図の例1．

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

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

図13-8: 窓関数による線形位相 FIRフィルタのインディシャル応答の例1．

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

図13-9: 窓関数による線形位相 FIRフィルタのゲイン線図の例2．

図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 は整数とする．

$$\begin{array}{}\text{(12-1)}& p\left(t\right)=f\left(t\right)\sum _{n=0}^{N-1}\delta \left(t-nT\right)=\sum _{n=0}^{N-1}f\left(nT\right)\delta \left(t-nT\right)\end{array}$$

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

$$\begin{array}{}\text{(12-2)}& \begin{array}{}\begin{array}{r}{c}_k\\ \\ \\ \end{array}\begin{array}{l}=\frac{1}{2\pi }{\int }_0^{2\pi }p\left(t\right)e^{-jkt}dt\\ =\frac{1}{2\pi }{\int }_0^{2\pi }\left\{\sum _{n=0}^{N-1}f\left(nT\right)\delta \left(t-nT\right)\right\}{e}^{-jkt}dt\\ =\frac{1}{2\pi }\sum _{n=0}^{N-1}f\left(nT\right)\left({\int }_{-\pi }^{\pi }\delta \left(t-nT\right)\mathrm{<mspace\; width="0.1em"></mspace>}{e}^{-jkt}dt\right)\\ =\frac{1}{2\pi }\sum _{n=0}^{N-1}f\left(nT\right)e^{-jknT}=\frac{1}{2\pi }\sum _{n=0}^{N-1}f\left(nT\right)e^{-j\frac{2\pi kn}{N}}\end{array}\end{array}\end{array}$$

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

$$\begin{array}{}\text{(12-3)}& F\left[k\right]=2\pi c_k=\sum _{n=0}^{N-1}f\left[n\right]e^{-j\frac{2\pi kn}{N}}\end{array}$$

以上から明らかな通り，離散フーリエ変換 (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) と呼ばれる．

$$\begin{array}{}\text{(12-4)}& \begin{array}{}\begin{array}{c}F\left[k\right]=\sum _{n=0}^{N-1}f\left[n\right]{W_N}^{kn},\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{<mspace\; width="0.1em"></mspace>}\left(0\le k\le N-1\right)\\ f\left[n\right]=\frac{1}{N}\sum _{k=0}^{N-1}F\left[k\right]{W_N}^{-kn},\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{<mspace\; width="0.1em"></mspace>}\left(0\le n\le N-1\right)\\ W_N=e^{-j\frac{2\pi }{N}}\end{array}\end{array}\end{array}$$

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

$$\begin{array}{}\text{(12-5)}& \begin{array}{c}\begin{array}{r}f\left[n\right]\\ \\ \\ \end{array}\begin{array}{l}=\frac{1}{N}\sum _{k=0}^{N-1}F\left[k\right]{W_N}^{-kn}=\frac{1}{N}\sum _{k=0}^{N-1}F\left[k\right]e^{j\frac{2\pi kn}{N}}\\ =\frac{F\left[0\right]}{N}+\frac{F\left[N/2\right]}{N}{e}^{jn\pi }+\frac{1}{N}\sum _{k=1}^{\left(N/2\right)-1}\left(F\left[k\right]e^{j\frac{2\pi kn}{N}}+\overline{F\left[k\right]}{e}^{-j\frac{2\pi kn}{N}}\right)\\ =\frac{F\left[0\right]}{N}+\frac{F\left[N/2\right]}{N}\mathrm{c}\mathrm{o}\mathrm{s}\left(n\pi \right)+\frac{1}{N}\sum _{k=1}^{\left(N/2\right)-1}\frac{j2F\left[k\right]e^{j\frac{2\pi kn}{N}}-\left(-j2\overline{F\left[k\right]}{e}^{-j\frac{2\pi kn}{N}}\right)}{j2}\\ =\frac{F\left[0\right]}{N}+\frac{F\left[N/2\right]}{N}\mathrm{s}\mathrm{i}\mathrm{n}\left(n\pi +\frac{\pi }{2}\right)+\sum _{k=1}^{\left(N/2\right)-1}{A}_k\frac{e^{j\left(\frac{2\pi kn}{N}+{\phi }_k\right)}-e^{-j\left(\frac{2\pi kn}{N}+{\phi }_k\right)}}{j2}\end{array}\\ A_k=\left|\frac{2F\left[k\right]}{N}\right|,\mathrm{<mspace\; width="0.4em"></mspace>}\mathrm{<mspace\; width="0.4em"></mspace>}{\phi }_k=\mathrm{A}\mathrm{r}\mathrm{g}\mathrm{<mspace\; width="0.1em"></mspace>}\left(j\frac{2F\left[k\right]}{N}\right),\mathrm{<mspace\; width="0.4em"></mspace>}\mathrm{<mspace\; width="0.4em"></mspace>}\left(1\le k\le \frac{N}{2}-1\right)\end{array}\end{array}$$

上記は (1) 角周波数を 2πk / N，振幅を A_k，初期位相を φ_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) の主要な定理．

線 形 性	$$\mathcal{F}\left\{a\mathrm{<mspace\; width="0.1em"></mspace>}f\left[n\right]+b\mathrm{<mspace\; width="0.1em"></mspace>}g\left[n\right]\right\}=a\mathrm{<mspace\; width="0.1em"></mspace>}\mathcal{F}\left\{f\left[n\right]\right\}+b\mathrm{<mspace\; width="0.1em"></mspace>}\mathcal{F}\left\{g\left[n\right]\right\}=a\mathrm{<mspace\; width="0.1em"></mspace>}F\left[k\right]+b\mathrm{<mspace\; width="0.1em"></mspace>}G\left[k\right]$$
畳み込み	$$\mathcal{F}\left\{\sum _{m=0}^{N-1}f\left[n-m\right]g\left[m\right]\right\}=\mathcal{F}\left\{f\left[n\right]\right\}\mathrm{<mspace\; width="0.1em"></mspace>}\mathcal{F}\left\{g\left[n\right]\right\}=F\left[k\right]G\left[k\right]$$
パーセバルの定理	$$\sum _{n=0}^{N-1}{\left|f\left[n\right]\right|}^2=\frac{1}{N}\sum _{k=0}^{N-1}{\left|F\left[k\right]\right|}^2$$
時間軸のシフト	$$\mathcal{F}\left\{f\left[n-a\right]\right\}=e^{-j\frac{2\pi ka}{N}}\mathcal{F}\left\{f\left[n\right]\right\}={W_N}^{ka}F\left[k\right]$$
周波数軸のシフト	$$\mathcal{F}\left\{e^{-j\frac{2\pi na}{N}}f\left[n\right]\right\}=\mathcal{F}\left\{{W_N}^{na}f\left[n\right]\right\}=F\left[k+a\right]$$

現実の信号分析と窓関数

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

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

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

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

$$\begin{array}{}\text{(12-6)}& w\left[n\right]=0.54-0.46\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{2\pi n}{N-1}\right),\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{<mspace\; width="0.1em"></mspace>}\left(0\le n\le N-1\right)\end{array}$$

上記のハミング窓を離散時間信号 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) は以下の式で記述される．

$$\begin{array}{}\text{(12-7)}& p\left(t\right)=f\left(t\right)\sum _{n=-\infty }^{\infty }\delta \left(t-nT\right)=\sum _{n=-\infty }^{\infty }f\left(nT\right)\delta \left(t-nT\right)\end{array}$$

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

$$\begin{array}{}\text{(12-8)}& \begin{array}{}\begin{array}{r}P\left(\omega \right)\\ \\ \\ \end{array}\begin{array}{l}={\int }_{-\infty }^{\infty }p\left(t\right)e^{-j\omega t}dt\\ ={\int }_{-\infty }^{\infty }\left\{\sum _{n=-\infty }^{\infty }f\left(nT\right)\delta \left(t-nT\right)\right\}{e}^{-j\omega t}dt\\ =\sum _{n=-\infty }^{\infty }f\left(nT\right)\left({\int }_{-\infty }^{\infty }\delta \left(t-nT\right){\mathrm{<mspace\; width="0.1em"></mspace>}e}^{-j\omega t}dt\right)\\ =\sum _{n=-\infty }^{\infty }f\left(nT\right)e^{-jn\omega T}\end{array}\end{array}\end{array}$$

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

$$\begin{array}{}\text{(12-9)}& P\left(\omega \right)=\sum _{n=-\infty }^{\infty }f\left[n\right]e^{-jn\omega T}\end{array}$$

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

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

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

$$\begin{array}{}\text{(12-10)}& F\left(\omega \right)=\sum _{n=-\infty }^{\infty }f\left[n\right]e^{-jn\omega T}\end{array}$$

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

$$\begin{array}{}\text{(12-11)}& f\left[n\right]=\frac{T}{2\pi }{\int }_{-\pi /T}^{\pi /T}F\left(\omega \right)\mathrm{<mspace\; width="0.1em"></mspace>}{e}^{jn\omega T}d\omega \end{array}$$

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

$$\begin{array}{}\text{(12-12)}& \begin{array}{c}f\left[n\right]=\frac{T}{2\pi }{\int }_{-\frac{\pi }{T}}^{\frac{\pi }{T}}F\left(\omega \right)\mathrm{<mspace\; width="0.1em"></mspace>}{e}^{jn\omega T}d\omega \\ \begin{array}{r}\therefore f\left[n\right]-\frac{TF\left(0\right)}{2\pi }\\ \\ \end{array}\begin{array}{l}~\frac{T}{2\pi }{\int }_{0+}^{\frac{\pi }{T}}F\left(\omega \right)\mathrm{<mspace\; width="0.1em"></mspace>}{e}^{jn\omega T}+\overline{F\left(\omega \right)}\mathrm{<mspace\; width="0.1em"></mspace>}{e}^{jn\omega T}d\omega \\ =\frac{T}{\pi }{\int }_{0+}^{\frac{\pi }{T}}\frac{jF\left(\omega \right)e^{jn\omega T}-\left\{-j\overline{F\left(\omega \right)}{e}^{jn\omega T}\right\}}{j2}d\omega \\ ={\int }_{0+}^{\frac{\pi }{T}}A\left(\omega \right)\frac{e^{j\left\{n\omega T+\phi \left(\omega \right)\right\}}-e^{-j\left\{n\omega T+\phi \left(\omega \right)\right\}}}{j2}d\omega \end{array}\\ A\left(\omega \right)=\left|\frac{T\mathrm{<mspace\; width="0.1em"></mspace>}F\left(\omega \right)}{\pi }\right|,\mathrm{<mspace\; width="0.4em"></mspace>}\mathrm{<mspace\; width="0.4em"></mspace>}\phi \left(\omega \right)=\mathrm{A}\mathrm{r}\mathrm{g}\mathrm{<mspace\; width="0.1em"></mspace>}\left\{j\frac{T\mathrm{<mspace\; width="0.1em"></mspace>}F\left(\omega \right)}{\pi }\right\},\mathrm{<mspace\; width="0.4em"></mspace>}\mathrm{<mspace\; width="0.4em"></mspace>}\left(0<\omega \le \frac{\pi }{T}\right)\end{array}\end{array}$$

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

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

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

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

$$\begin{array}{}\text{(12-13)}& F\left(\omega \right)=\frac{1}{T}\sum _{n=-\infty }^{\infty }{F}_a\left(\omega -\frac{2\pi n}{T}\right)\end{array}$$

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

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

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

表12-2: 離散時間フーリエ変換 (DTFT) の主要な定理．

線 形 性	$$\mathcal{F}\left\{a\mathrm{<mspace\; width="0.1em"></mspace>}f\left[n\right]+b\mathrm{<mspace\; width="0.1em"></mspace>}g\left[n\right]\right\}=a\mathrm{<mspace\; width="0.1em"></mspace>}\mathcal{F}\left\{f\left[n\right]\right\}+b\mathrm{<mspace\; width="0.1em"></mspace>}\mathcal{F}\left\{g\left[n\right]\right\}=a\mathrm{<mspace\; width="0.1em"></mspace>}F\left(\omega \right)+b\mathrm{<mspace\; width="0.1em"></mspace>}G\left(\omega \right)$$
畳み込み	$$\mathcal{F}\left\{\sum _{k=-\infty }^{\infty }f\left[n-k\right]g\left[k\right]\right\}=\mathcal{F}\left\{f\left[n\right]\right\}\mathrm{<mspace\; width="0.1em"></mspace>}\mathcal{F}\left\{g\left[n\right]\right\}=F\left(\omega \right)G\left(\omega \right)$$
パーセバルの定理	$$\sum _{n=-\infty }^{\infty }{\left|f\left[n\right]\right|}^2=\frac{1}{2\pi }{\int }_{-\pi /T}^{\pi /T}{\left|F\left(\omega \right)\right|}^{2}d\omega$$
時間軸のシフト	$$\mathcal{F}\left\{f\left[n-a\right]\right\}=e^{-ja\omega T}\mathcal{F}\left\{f\left[n\right]\right\}=e^{-ja\omega T}F\left(\omega \right)$$
周波数軸のシフト	$$\mathcal{F}\left\{e^{-jat}f\left[n\right]\right\}=F\left(\omega +a\right)$$

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

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

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

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

$$\begin{array}{}\text{(12-14)}& \mathcal{F}\left[{\int }_{-\infty }^{\infty }f\left(t-\tau \right)g\left(\tau \right)d\tau \right]=F\left(\omega \right)G\left(\omega \right)\end{array}$$

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

$$\begin{array}{}\text{(12-15)}& \mathcal{F}\left[\frac{1}{2\pi }{\int }_{-\infty }^{\infty }F\left(\omega -\tau \right)G\left(\tau \right)d\tau \right]=f\left(t\right)g\left(t\right)\end{array}$$

櫛型関数のフーリエ級数

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

$$\begin{array}{}\text{(12-16)}& \begin{array}{}\begin{array}{r}{c}_n\\ \\ \\ \end{array}\begin{array}{l}=\frac{1}{T}{\int }_{-T/2}^{T/2}\left(\sum _{n=-\infty }^{\infty }\delta \left(t-nT\right)\mathrm{<mspace\; width="0.1em"></mspace>}{e}^{-j\frac{2\pi nt}{T}}\right)dt\\ =\frac{1}{T}\sum _{n=-\infty }^{\infty }\left({\int }_{-T/2}^{T/2}\delta \left(t-nT\right)\mathrm{<mspace\; width="0.1em"></mspace>}{e}^{-j\frac{2\pi nt}{T}}dt\right)\\ =\frac{1}{T}\sum _{n=-\infty }^{\infty }\left({\int }_{-T/2}^{T/2}\delta \left(t\right)\mathrm{<mspace\; width="0.1em"></mspace>}{e}^{-j\frac{2\pi nt}{T}}dt\right)\\ =\frac{1}{T}\sum _{n=-\infty }^{\infty }{e}^0=\frac{1}{T}\end{array}\end{array}\end{array}$$

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

$$\begin{array}{}\text{(12-17)}& \sum _{n=-\infty }^{\infty }\delta \left(t-nT\right)=\frac{1}{T}\sum _{n=-\infty }^{\infty }{e}^{j\frac{2\pi nt}{T}}\end{array}$$

櫛型関数のフーリエ変換

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

$$\begin{array}{}\text{(12-18)}& \mathcal{F}\left\{\sum _{n=-\infty }^{\infty }\delta \left(t-nT\right)\right\}=\sum _{n=-\infty }^{\infty }{e}^{-jn\omega T}\end{array}$$

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

$$\begin{array}{}\text{(12-19)}& \begin{array}{}\begin{array}{r}\mathcal{F}\left\{\sum _{n=-\infty }^{\infty }\delta \left(t-nT\right)\right\}\\ \\ \\ \\ \end{array}\begin{array}{l}={\int }_{-\infty }^{\infty }\left\{\sum _{n=-\infty }^{\infty }\delta \left(t-nT\right)\right\}{e}^{-j\omega t}dt\\ =\frac{1}{T}{\int }_{-\infty }^{\infty }\left(\sum _{n=-\infty }^{\infty }{e}^{j\frac{2\pi nt}{T}}\right)e^{-j\omega t}dt\\ =\frac{1}{T}\sum _{n=-\infty }^{\infty }\left({\int }_{-\infty }^{\infty }{e}^{j\frac{2\pi nt}{T}}{e}^{-j\omega t}dt\right)\\ =\frac{1}{T}\sum _{n=-\infty }^{\infty }\left({\int }_{-\infty }^{\infty }{e}^{-j\left(\omega -\frac{2\pi n}{T}\right)t}dt\right)\\ =\frac{2\pi }{T}\sum _{n=-\infty }^{\infty }\delta \left(\omega -\frac{2\pi n}{T}\right)\end{array}\end{array}\end{array}$$

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

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

$$\begin{array}{}\text{(12-20)}& p\left(t\right)=f\left(t\right)\sum _{n=-\infty }^{\infty }\delta \left(t-nT\right)\end{array}$$

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

$$\begin{array}{}\text{(12-21)}& \begin{array}{}\begin{array}{r}F\left(\omega \right)\\ \\ \end{array}\begin{array}{l}=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }{F}_a\left(\tau \right)\frac{2\pi }{T}\left\{\sum _{n=-\infty }^{\infty }\delta \left(\omega -\frac{2\pi n}{T}-\tau \right)\right\}d\tau \\ =\frac{1}{T}\sum _{n=-\infty }^{\infty }\left\{{\int }_{-\infty }^{\infty }{F}_a\left(\tau \right)\delta \left(\omega -\frac{2\pi n}{T}-\tau \right)d\tau \right\}\\ =\frac{1}{T}\sum _{n=-\infty }^{\infty }{F}_a\left(\omega -\frac{2\pi n}{T}\right)\end{array}\end{array}\end{array}$$

本章の参考文献

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

11. IIR フィルタの概要

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

11.1. IIR フィルタの構造

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

$$\begin{array}{}\text{(11-1)}& H\left(z\right)=\frac{Y\left(z\right)}{X\left(z\right)}=\frac{{\sum }_{r=0}^M{b}_r\mathrm{<mspace\; width="0.1em"></mspace>}{z}^{-r}}{1-\sum _{k=1}^N{a}_k{z}^{-k}}\end{array}$$

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

$$\begin{array}{}\text{(11-2)}& y\left[n\right]=\sum _{k=1}^N{a}_k\mathrm{<mspace\; width="0.1em"></mspace>}y\left[n-k\right]+\sum _{r=0}^M{b}_r\mathrm{<mspace\; width="0.1em"></mspace>}x\left[n-r\right]\end{array}$$

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

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

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

図11-1: IIR フィルタの構造 (直接形)．

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

$$\begin{array}{}\text{(11-3)}& H\left(z\right)=\frac{{\sum }_{r=0}^M{b}_r\mathrm{<mspace\; width="0.1em"></mspace>}{z}^{-r}}{1-\sum _{k=1}^N{a}_k{z}^{-k}}=A\prod _{k=1}^{⌊\left(N+1\right)/2⌋}\frac{1+{\beta }_{1k}{z}^{-1}+{\beta }_{2k}{z}^{-2}}{1-{\alpha }_{1k}{z}^{-1}-{\alpha }_{2k}{z}^{-2}},\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{<mspace\; width="0.1em"></mspace>}\left(M\le N\right)\end{array}$$

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

図11-2: IIR フィルタの構造 (縦続形)．

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

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

$$\begin{array}{}\text{(11-4)}& h\mathrm{<mspace\; width="0.1em"></mspace>}\left[n\right]=g\mathrm{<mspace\; width="0.1em"></mspace>}\left(nT\right),\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{<mspace\; width="0.1em"></mspace>}\left(n\ge 0\right)\end{array}$$

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

システム伝達関数

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

$$\begin{array}{}\text{(11-5)}& G\left(s\right)=\frac{1}{\tau s+1}=\frac{{\omega }_c}{s+{\omega }_c}\end{array}$$

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

$$\begin{array}{}\text{(11-6)}& g\left(t\right)=\frac{1}{\tau }{e}^{-\frac{t}{\tau }},\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{<mspace\; width="0.1em"></mspace>}\left(t\ge 0\right)\end{array}$$

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

$$\begin{array}{}\text{(11-7)}& H\left(z\right)=\frac{{\omega }_c}{1-e^{-{\omega }_{c}T}{\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}}\end{array}$$

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

$$\begin{array}{}\text{(11-8)}& h\left[n\right]=\frac{1}{\tau }{e}^{-\frac{nT}{\tau }},\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{<mspace\; width="0.1em"></mspace>}\left(n\ge 0\right)\end{array}$$

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

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

$$\begin{array}{}\text{(11-9)}& H_N\left(z\right)=\frac{\underset{s\to 0}{\lim }G\left(s\right)}{\underset{z\to 1}{\lim }H\left(z\right)}H\left(z\right)=\frac{1-e^{-{\omega }_{c}T}}{1-e^{-{\omega }_{c}T}{\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}}\end{array}$$

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

$$\begin{array}{}\text{(11-10)}& H_N\left(e^{j\omega T}\right)=\frac{1-e^{-{\omega }_{c}T}}{1-e^{-{\omega }_{c}T}{e}^{-j\omega T}}\end{array}$$

差分方程式

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

$$\begin{array}{}\text{(11-11)}& \begin{array}{}\begin{array}{c}{H}_N\left(z\right)=\frac{Y\left(z\right)}{X\left(z\right)}=\frac{1-e^{-{\omega }_{c}T}}{1-e^{-{\omega }_{c}T}{\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}}\\ \therefore \left(1-e^{-{\omega }_{c}T}{\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}\right)\mathrm{<mspace\; width="0.1em"></mspace>}Y\left(z\right)=\left(1-e^{-{\omega }_{c}T}\right)\mathrm{<mspace\; width="0.1em"></mspace>}X\left(z\right)\\ y\left[n\right]-e^{-{\omega }_{c}T}\mathrm{<mspace\; width="0.1em"></mspace>}y\left[n-1\right]=\left(1-e^{-{\omega }_{c}T}\right)\mathrm{<mspace\; width="0.1em"></mspace>}x\left[n\right]\\ \therefore y\left[n\right]=e^{-{\omega }_{c}T}\mathrm{<mspace\; width="0.1em"></mspace>}y\left[n-1\right]+\left(1-e^{-{\omega }_{c}T}\right)\mathrm{<mspace\; width="0.1em"></mspace>}x\left[n\right]\end{array}\end{array}\end{array}$$

ボード線図

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

図11-3: インパルス不変設計法によるローパスフィルタのボード線図の例．

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

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

図11-3: インパルス不変設計法によるローパスフィルタのボード線図の例．

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

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

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

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

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

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

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

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

$$\begin{array}{}\text{(11-12)}& {\omega }_a=\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{\omega T}{2}\right)\end{array}$$

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

図11-6: 双一次変換による設計法の原理．

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

$$\begin{array}{}\text{(11-13)}& j{\omega }_a=j\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{\omega T}{2}\right)=j\frac{\sin \left(\frac{\omega T}{2}\right)}{\cos \left(\frac{\omega T}{2}\right)}=j\frac{\frac{e^{j\frac{\omega T}{2}}-e^{-j\frac{\omega T}{2}}}{2j}}{\frac{e^{j\frac{\omega T}{2}}+e^{-j\frac{\omega T}{2}}}{2}}=\frac{e^{j\omega T}-1}{e^{j\omega T}+1}\end{array}$$

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

$$\begin{array}{}\text{(11-14)}& s_a=\frac{e^{sT}-1}{e^{sT}+1}\end{array}$$

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

$$\begin{array}{}\text{(11-15)}& s_a=\frac{e^{sT}-1}{e^{sT}+1}=\frac{z-1}{z+1}=\frac{1-{\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}}{1+{\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}}\end{array}$$

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

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

$$\begin{array}{}\text{(11-16)}& G\left(s\right)=\frac{1}{\tau s+1}=\frac{{\omega }_c}{s+{\omega }_c}\end{array}$$

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

$$\begin{array}{}\text{(11-17)}& {\omega }_a=\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{{\omega }_{c}T}{2}\right)\end{array}$$

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

$$\begin{array}{}\text{(11-18)}& H\left(z\right)=\frac{{\omega }_a}{s+{\omega }_a}=\frac{{\omega }_a}{\frac{1-{\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}}{1+{\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}}+{\omega }_a}=\frac{{\omega }_a\left(1+{\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}\right)}{\left(1+{\omega }_a\right)-\left(1-{\omega }_a\right){\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}}\end{array}$$

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

$$\begin{array}{}\text{(11-19)}& H\left(e^{j\omega T}\right)=\frac{{\omega }_a\left(1+e^{-j\omega T}\right)}{\left(1+{\omega }_a\right)-\left(1-{\omega }_a\right)e^{-j\omega T}}\end{array}$$

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

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

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

$$\begin{array}{}\text{(11-20)}& \begin{array}{}\begin{array}{c}H\left(z\right)=\frac{Y\left(z\right)}{X\left(z\right)}=\frac{{\omega }_a\left(1+{\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}\right)}{\left(1+{\omega }_a\right)-\left(1-{\omega }_a\right){\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}}\\ \therefore \left\{\left(1+{\omega }_a\right)-\left(1-{\omega }_a\right){\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}\right\}Y\left(z\right)={\omega }_a\left(1+{\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}\right)X\left(z\right)\\ \left(1+{\omega }_a\right)\mathrm{<mspace\; width="0.1em"></mspace>}y\left[n\right]-\left(1-{\omega }_a\right)\mathrm{<mspace\; width="0.1em"></mspace>}y\left[n-1\right]={\omega }_a\left(x\left[n\right]+x\left[n-1\right]\right)\\ \therefore y\left[n\right]=\frac{1-{\omega }_a}{1+{\omega }_a}\mathrm{<mspace\; width="0.1em"></mspace>}y\left[n-1\right]+\frac{{\omega }_a}{1+{\omega }_a}\left(x\left[n\right]+x\left[n-1\right]\right)\end{array}\end{array}\end{array}$$

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

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

図11-7: 双一次変換による設計法におけるローパスフィルタのボード線図の例．

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

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

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

$$\begin{array}{}\text{(11-21)}& G\left(s\right)=\frac{\tau s}{\tau s+1}=1-\frac{1}{\tau s+1}=1-\frac{{\omega }_c}{s+{\omega }_c}\end{array}$$

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

$$\begin{array}{}\text{(11-22)}& {\omega }_a=\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{{\omega }_{c}T}{2}\right)\end{array}$$

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

$$\begin{array}{}\text{(11-23)}& H\left(z\right)=1-\frac{{\omega }_a}{s+{\omega }_a}=1-\frac{{\omega }_a}{\frac{1-{\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}}{1+{\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}}+{\omega }_a}=1-\frac{{\omega }_a\left(1+{\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}\right)}{\left(1+{\omega }_a\right)-\left(1-{\omega }_a\right){\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}}\end{array}$$

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

$$\begin{array}{}\text{(11-24)}& H\left(e^{j\omega T}\right)=1-\frac{{\omega }_a\left(1+e^{-j\omega T}\right)}{\left(1+{\omega }_a\right)-\left(1-{\omega }_a\right)e^{-j\omega T}}\end{array}$$

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

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

$$\begin{array}{}\text{(11-25)}& \begin{array}{}\begin{array}{c}H\left(z\right)=\frac{Y\left(z\right)}{X\left(z\right)}=1-\frac{{\omega }_a\left(1+{\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}\right)}{\left(1+{\omega }_a\right)-\left(1-{\omega }_a\right){\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}}\\ \therefore \left\{\left(1+{\omega }_a\right)-\left(1-{\omega }_a\right){\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}\right\}Y\left(z\right)=\left(1-{\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}\right)X\left(z\right)\\ \left(1+{\omega }_a\right)\mathrm{<mspace\; width="0.1em"></mspace>}y\left[n\right]-\left(1-{\omega }_a\right)\mathrm{<mspace\; width="0.1em"></mspace>}y\left[n-1\right]=x\left[n\right]-x\left[n-1\right]\\ \therefore y\left[n\right]=\frac{1-{\omega }_a}{1+{\omega }_a}\mathrm{<mspace\; width="0.1em"></mspace>}y\left[n-1\right]+\frac{1}{1+{\omega }_a}\left(x\left[n\right]-x\left[n-1\right]\right)\end{array}\end{array}\end{array}$$

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

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

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

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

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

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

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

$$\begin{array}{}\text{(11-26)}& G\left(s\right)=\frac{\tau s}{\tau s+1}=1-\frac{1}{\tau s+1}\end{array}$$

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

$$\begin{array}{}\text{(11-27)}& g\left(t\right)=\delta \left(t\right)-\frac{1}{\tau }{e}^{-\frac{t}{\tau }},\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{<mspace\; width="0.1em"></mspace>}\left(t\ge 0\right)\end{array}$$

上記から明らかな通り，ハイパスフィルタ及びノッチフィルタにおける伝達関数の分子と分母の次数は等しいため，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 に示す．上のブロックダイヤグラムは論理モデルを，下のブロックダイヤグラムは実装モデルを表す．

図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 変換器の出力が階段関数となる事から，階段関数の信号における各階段の高さを示す係数を意味する．

図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 とする．

図10-3: 長方形とインパルス．

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

$$\begin{array}{}\text{(10-1)}& x\left[n\right]{\int }_{-\infty }^{\infty }\delta \left(t-nT\right)dt=x\left[n\right]\end{array}$$

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

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

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

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

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

$$\begin{array}{}\text{(10-2)}& TH\left(z\right)=\frac{{\omega }_{c}T}{1-e^{-{\omega }_{c}T}{\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}}\end{array}$$

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

$$\begin{array}{}\text{(10-3)}& TH\left(e^{j\omega T}\right)=\frac{{\omega }_{c}T}{1-e^{-{\omega }_{c}T}{e}^{-j\omega T}}\end{array}$$

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

図10-4: ゲインを補正ローパスフィルタのゲイン線図の例．

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

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

$$\begin{array}{}\text{(10-4)}& \underset{\omega \to 0}{\lim }T\mathrm{<mspace\; width="0.1em"></mspace>}H\left(e^{j\omega T}\right)=\underset{\omega \to 0}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{{\omega }_{c}T}{1-e^{-{\omega }_{c}T}{e}^{-j\omega T}}=\frac{{\omega }_{c}T}{1-e^{-{\omega }_{c}T}}\end{array}$$

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

図10-5: 階段関数による連続時間信号の近似．

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

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

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

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

$$\begin{array}{}\text{(10-5)}& {\int }_0^{\infty }g\left(t\right)\mathit{dt}=<mspace\; width="0.1em"></mspace>\underset{s\to 0}{\lim }G\left(s\right)\end{array}$$

上記は一見不思議な式であるが，式 (5-1) に示したラプラス変換の定義から以下の通り求められる．

$$\begin{array}{}\text{(10-6)}& \underset{s\to 0}{\lim }F\left(s\right)=\underset{s\to 0}{\lim }{\int }_0^{\infty }f\left(t\right)e^{-st}dt={\int }_0^{\infty }f\left(t\right)dt\end{array}$$

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

$$\begin{array}{}\text{(10-7)}& \sum _{n=0}^{\infty }T\mathrm{<mspace\; width="0.1em"></mspace>}h\left[n\right]=T\mathrm{<mspace\; width="0.1em"></mspace>}\underset{z\to 1}{\lim }H\left(z\right)\end{array}$$

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

$$\begin{array}{}\text{(10-8)}& H_N\left(z\right)=\frac{\underset{s\to 0}{\lim }G\left(s\right)}{T\underset{z\to 1}{\lim }H\left(z\right)}TH\left(z\right)=\frac{\underset{s\to 0}{\lim }G\left(s\right)}{\underset{z\to 1}{\lim }H\left(z\right)}H\left(z\right)\end{array}$$

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

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

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

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

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

$$\begin{array}{}\text{(10-9)}& G\left(s\right)=\frac{1}{\tau s+1}=\frac{{\omega }_c}{s+{\omega }_c}\end{array}$$

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

$$\begin{array}{}\text{(10-10)}& H\left(z\right)=\frac{{\omega }_c}{1-e^{-{\omega }_{c}T}{\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}}\end{array}$$

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

$$\begin{array}{}\text{(10-11)}& H_N\left(z\right)=\frac{\underset{s\to 0}{\lim }G\left(s\right)}{\underset{z\to 1}{\lim }H\left(z\right)}H\left(z\right)=\frac{1-e^{-{\omega }_{c}T}}{1-e^{-{\omega }_{c}T}{\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}}\end{array}$$

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

$$\begin{array}{}\text{(10-12)}& H_N\left(e^{j\omega T}\right)=\frac{1-e^{-{\omega }_{c}T}}{1-e^{-{\omega }_{c}T}{e}^{-j\omega T}}\end{array}$$

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

図10-6: 正規化されたローパスフィルタのゲイン線図の例．

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

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

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

$$\begin{array}{}\text{(10-13)}& G\left(s\right)=\frac{1}{s}-\frac{e^{-sT}}{s}=\frac{1-e^{-sT}}{s}\end{array}$$

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

$$\begin{array}{}\text{(10-14)}& G_N\left(s\right)=\frac{1}{T}G\left(s\right)=\frac{1-e^{-sT}}{sT}\end{array}$$

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

$$\begin{array}{}\text{(10-15)}& G_N\left(j\omega \right)=\frac{1-e^{-j\omega T}}{j\omega T}\end{array}$$

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

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

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

図10-7: 0 次ホールドのボード線図．

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

本章の参考文献

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

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

2017/5/29 10.4 を一部修正．

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

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

9.1. サンプリング定理

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

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

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

9.2. サンプル値の表現

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

図9-1: サンプラーの論理モデル．

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

$$\begin{array}{}\text{(9-1)}& \begin{array}{}\begin{array}{r}p\left(t\right)\\ \\ \end{array}\begin{array}{l}=x\left(t\right)\mathrm{<mspace\; width="0.1em"></mspace>}\left\{\delta \left(t\right)+\delta \left(t-T\right)+\delta \left(t-2T\right)+\delta \left(t-3T\right)+\mathrm{<mspace\; width="0.1em"></mspace>}\cdots \right\}\\ =x\left(t\right)\sum _{n=0}^{\infty }\delta \left(t-nT\right)\\ =\sum _{n=0}^{\infty }x\left(nT\right)\delta \left(t-nT\right),\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{<mspace\; width="0.1em"></mspace>}\left(t\ge 0\right)\end{array}\end{array}\end{array}$$

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

9.3. Z 変換の導出

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

$$\begin{array}{}\text{(9-2)}& P\left(s\right)=\sum _{n=0}^{\infty }x\left(nT\right)e^{-nTs}\end{array}$$

上記の式において z = e^Ts に置き換えると，以下の通り Z 変換が得られる．

$$\begin{array}{}\text{(9-3)}& P\left(z\right)=\sum _{n=0}^{\infty }x\left(nT\right)z^{-n}\end{array}$$

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

$$\begin{array}{}\text{(9-4)}& \begin{array}{}\begin{array}{c}p\left(t\right)=\sum _{n=0}^{\infty }x\left[n\right]\delta \left(t-nT\right),\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{<mspace\; width="0.1em"></mspace>}\left(t\ge 0\right)\\ P\left(z\right)=\sum _{n=0}^{\infty }x\left[n\right]z^{-n}\end{array}\end{array}\end{array}$$

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

9.4. Z 変換の定義

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

$$\begin{array}{}\text{(9-5)}& F\left(z\right)=\sum _{n=0}^{\infty }f\left[n\right]z^{-n}\end{array}$$

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

$$\begin{array}{}\text{(9-6)}& F\left(z\right)=\sum _{n=0}^{\infty }f\left[n\right]z^{-n}=3+z^{-1}+4z^{-2}\end{array}$$

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

$$\begin{array}{}\text{(9-7)}& Z\left\{f\left[n\right]\right\}=F\left(z\right)\end{array}$$

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

$$\begin{array}{}\text{(9-8)}& Z^{-1}\left\{F\left(z\right)\right\}=f\left[n\right]\end{array}$$

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

$$\begin{array}{}\text{(9-9)}& f\left[n\right]=\frac{1}{j2\pi }{\oint }_C^{}F\left(z\right)\mathrm{<mspace\; width="0.1em"></mspace>}{z}^{n-1}dz\end{array}$$

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

9.5. Z 変換の公式と定理

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

表9-1: Z変換の主要な公式．

項番	$f\left[n\right],\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{<mspace\; width="0.1em"></mspace>}\left(n\ge 0\right)$	$F\left(z\right)$	収束域
1	$$\delta \left[n\right]=\left\{\begin{array}{c}1,\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{<mspace\; width="0.1em"></mspace>}\left(n=0\right)\\ 0,\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{<mspace\; width="0.1em"></mspace>}\left(n\ne 0\right)\end{array}\right.$$	$$1$$
2	$$u\left[n\right]=\left\{\begin{array}{c}1,\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{<mspace\; width="0.1em"></mspace>}\left(n\ge 0\right)\\ 0,\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{<mspace\; width="0.1em"></mspace>}\left(n<0\right)\end{array}\right.$$	$$\frac{1}{{1-z}^{-1}}$$	$$\left|z\right|>1$$
3	$$1$$	$$\frac{1}{{1-z}^{-1}}$$	$$\left|z\right|>1$$
4	$$n$$	$$\frac{z^{-1}}{{\left({1-z}^{-1}\right)}^2}$$	$$\left|z\right|>1$$
5	$${\alpha }^n$$	$$\frac{1}{{1-\alpha \mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}}$$	$$\left|z\right|>\left|\alpha \right|$$
6	$$n{\mathrm{<mspace\; width="0.1em"></mspace>}\alpha }^n$$	$$\frac{{\alpha \mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}}{{\left({1-\alpha \mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}\right)}^2}$$	$$\left|z\right|>\left|\alpha \right|$$

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

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

$$\begin{array}{}\text{(9-10)}& \sum _{n=0}^{\infty }a{r}^n=\frac{a}{1-r},\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{<mspace\; width="0.1em"></mspace>}\left|r\right|<1\end{array}$$

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

表9-2: Z変換の主要な定理．

線 形 性	$$Z\left\{a\mathrm{<mspace\; width="0.1em"></mspace>}f\left[n\right]+b\mathrm{<mspace\; width="0.1em"></mspace>}g\left[n\right]\right\}=a\mathrm{<mspace\; width="0.1em"></mspace>}Z\left\{f\left[n\right]\right\}+b\mathrm{<mspace\; width="0.1em"></mspace>}Z\left\{g\left[n\right]\right\}=a\mathrm{<mspace\; width="0.1em"></mspace>}F\left(z\right)+b\mathrm{<mspace\; width="0.1em"></mspace>}G\left(z\right)$$
畳み込み	$$Z\left\{\sum _{k=0}^{\infty }f\left[n-k\right]\mathrm{<mspace\; width="0.1em"></mspace>}g\left[k\right]\right\}=Z\left\{f\left[n\right]\right\}\mathrm{<mspace\; width="0.1em"></mspace>}Z\left\{g\left[n\right]\right\}=F\left(z\right)G\left(z\right)$$
時間軸のシフト	$$Z\left\{f\left[n-k\right]\mathrm{<mspace\; width="0.1em"></mspace>}\right\}=z^{-k}\mathrm{<mspace\; width="0.1em"></mspace>}Z\left\{f\left[n\right]\right\}=z^{-k}\mathrm{<mspace\; width="0.1em"></mspace>}F\left(z\right),\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{<mspace\; width="0.1em"></mspace>}\left(k\ge 0\right)$$
指数数列の乗算	$$Z\left\{a^{n}f\left[n\right]\right\}=F\left(a^{-1}\mathrm{<mspace\; width="0.1em"></mspace>}z\right)$$
F (z) の微分	$$Z\left\{n\mathrm{<mspace\; width="0.1em"></mspace>}f\left[n\right]\right\}=-z\frac{d}{dz}F\left(z\right)$$
最終値定理	$$f\left[\mathrm{\infty }\right]<mspace\; width="0.1em"></mspace>=\underset{z\to 1}{\lim }\left({1-z}^{-1}\right)F\left(z\right)$$
初期値定理	$$f\left[0\right]<mspace\; width="0.1em"></mspace>=\underset{z\to \mathrm{\infty }}{\lim }F\left(z\right)$$

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

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

表9 3: ラプラス変換と Z 変換の対応関係．

原関数 $f\left(t\right),\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{<mspace\; width="0.1em"></mspace>}\left(t\ge 0\right)$	像関数 $F\left(s\right)$	$f\left[n\right],\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{<mspace\; width="0.1em"></mspace>}\left(n\ge 0\right)$	$F\left(z\right)$
$$u\left(t\right)=\left\{\begin{array}{c}1,\mathrm{<mspace\; width="0.1em"></mspace>}\left(t\ge 0\right)\\ 0,\mathrm{<mspace\; width="0.1em"></mspace>}\left(t<0\right)\end{array}\right.$$	$$\frac{1}{s}$$	$$u\left[n\right]=\left\{\begin{array}{c}1,\mathrm{<mspace\; width="0.1em"></mspace>}\left(n\ge 0\right)\\ 0,\mathrm{<mspace\; width="0.1em"></mspace>}\left(n<0\right)\end{array}\right.$$	$$\frac{1}{{1-z}^{-1}}$$
$$t$$	$$\frac{1}{s^2}$$	$$nT$$	$$\frac{T{z}^{-1}}{{\left({1-z}^{-1}\right)}^2}$$
$$e^{-at}$$	$$\frac{1}{s+a}$$	$$e^{-naT}$$	$$\frac{1}{{1-e^{-aT}z}^{-1}}$$

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

$$\begin{array}{}\text{(9-11)}& \mathfrak{L}\left[e^{-at}\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{c}\mathrm{o}\mathrm{s}\left(\omega t\right)\right]=\frac{s+a}{{\left(s+a\right)}^2+{\omega }^2}\end{array}$$

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

$$\begin{array}{}\text{(9-12)}& \frac{s+a}{{\left(s+a\right)}^2+{\omega }^2}=\frac{1}{2}\left\{\frac{1}{s+\left(a+j\omega \right)}+\frac{1}{s+\left(a-j\omega \right)}\right\}\end{array}$$

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

$$\begin{array}{}\text{(9-13)}& \begin{array}{}\begin{array}{r}Z\left\{\frac{s+a}{{\left(s+a\right)}^2+{\omega }^2}\right\}\\ \\ \\ \end{array}\begin{array}{l}=\frac{1}{2}\left(\frac{1}{{1-e^{-\left(a+j\omega \right)T}z}^{-1}}+\frac{1}{{1-e^{-\left(a-j\omega \right)T}z}^{-1}}\right)\\ =\frac{1}{2}\mathrm{<mspace\; width="0.1em"></mspace>}\frac{2-\left(e^{-\left(a+j\omega \right)T}+e^{-\left(a-j\omega \right)T}\right)\mathrm{<mspace\; width="0.1em"></mspace>}{z}^{-1}}{\left({1-e^{-\left(a+j\omega \right)T}z}^{-1}\right)\left({1-e^{-\left(a-j\omega \right)T}z}^{-1}\right)}\\ =\frac{1}{2}\frac{2-e^{-aT}\left(e^{-j\omega T}+e^{j\omega T}\right)z^{-1}}{1-e^{-aT}\left(e^{-j\omega T}+e^{j\omega T}\right)\mathrm{<mspace\; width="0.1em"></mspace>}{z}^{-1}+e^{-2aT}{z}^{-2}}\\ =\frac{1-e^{-aT}\mathrm{c}\mathrm{o}\mathrm{s}\left(\omega T\right){\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}}{1-2e^{-aT}{\mathrm{c}\mathrm{o}\mathrm{s}\left(\omega T\right)\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}+e^{-2aT}{z}^{-2}}\end{array}\end{array}\end{array}$$

別の例として，連続時間関数 e^−at sin (ωt) のラプラス変換を以下に示す．

$$\begin{array}{}\text{(9-14)}& \mathfrak{L}\left[e^{-at}\mathrm{<mspace\; width="0.1em"></mspace>}\mathrm{s}\mathrm{i}\mathrm{n}\left(\omega t\right)\right]=\frac{\omega }{{\left(s+a\right)}^2+{\omega }^2}\end{array}$$

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

$$\begin{array}{}\text{(9-15)}& Z\left\{\frac{\omega }{{\left(s+a\right)}^2+{\omega }^2}\right\}=\frac{e^{-aT}\mathrm{s}\mathrm{i}\mathrm{n}\left(\omega T\right){\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}}{1-2e^{-aT}{\mathrm{c}\mathrm{o}\mathrm{s}\left(\omega T\right)\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}+e^{-2aT}{z}^{-2}}\end{array}$$

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

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

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

$$\begin{array}{}\text{(9-16)}& y\mathrm{<mspace\; width="0.1em"></mspace>}\left[n\right]-a\mathrm{<mspace\; width="0.1em"></mspace>}y\mathrm{<mspace\; width="0.1em"></mspace>}\left[n-1\right]=x\mathrm{<mspace\; width="0.1em"></mspace>}\left[n\right]\end{array}$$

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

$$\begin{array}{}\text{(9-17)}& Y\left(z\right)-a\mathrm{<mspace\; width="0.1em"></mspace>}{z}^{-1}Y\left(z\right)=X\left(z\right)\end{array}$$

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

$$\begin{array}{}\text{(9-18)}& Y\left(z\right)=\frac{X\left(z\right)}{1-a{\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}}\end{array}$$

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

$$\begin{array}{}\text{(9-19)}& y\mathrm{<mspace\; width="0.1em"></mspace>}\left[n\right]=Z^{-1}\left\{\frac{X\left(z\right)}{1-a{\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}}\right\}\end{array}$$

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 の場合の時間軸シフトの定理を以下に示す．

$$\begin{array}{}\text{(9-20)}& Z\left\{f\left[n+2\right]\mathrm{<mspace\; width="0.1em"></mspace>}\right\}=z^2\mathrm{<mspace\; width="0.1em"></mspace>}Z\left\{f\left[n\right]\right\}-z^2\mathrm{<mspace\; width="0.1em"></mspace>}f\left[0\right]-z\mathrm{<mspace\; width="0.1em"></mspace>}f\left[1\right]=z^2\mathrm{<mspace\; width="0.1em"></mspace>}F\left(z\right)-z^2\mathrm{<mspace\; width="0.1em"></mspace>}f\left[0\right]-z\mathrm{<mspace\; width="0.1em"></mspace>}f\left[1\right]\end{array}$$

9.8. 逆 Z 変換

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

$$\begin{array}{}\text{(9-21)}& Y\left(z\right)=\frac{b_0+b_1\mathrm{<mspace\; width="0.1em"></mspace>}{z}^{-1}}{{a_0+a_1\mathrm{<mspace\; width="0.1em"></mspace>}{z}^{-1}+a}_2\mathrm{<mspace\; width="0.1em"></mspace>}{z}^{-2}}\end{array}$$

べき級数展開による方法

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

$$\begin{array}{}\text{(9-22)}& \begin{array}{}\begin{array}{r}Y\left(z\right)\\ \end{array}\begin{array}{l}=\frac{b_0+b_1\mathrm{<mspace\; width="0.1em"></mspace>}{z}^{-1}}{{a_0+a_1\mathrm{<mspace\; width="0.1em"></mspace>}{z}^{-1}+a}_2\mathrm{<mspace\; width="0.1em"></mspace>}{z}^{-2}}\\ =y\left[0\right]+y\left[1\right]\mathrm{<mspace\; width="0.1em"></mspace>}{z}^{-1}+y\left[2\right]\mathrm{<mspace\; width="0.1em"></mspace>}{z}^{-2}+\cdots \end{array}\end{array}\end{array}$$

部分分数分解による方法

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

$$\begin{array}{}\text{(9-23)}& Y\left(z\right)=k_0+\frac{k_1}{1-\alpha \mathrm{<mspace\; width="0.1em"></mspace>}{z}^{-1}}+\frac{k_2}{1-\beta {\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}}\end{array}$$

ここで，k₁, k₂ は以下の様に留数定理 (ヘビサイドの展開定理) から，k₀ は係数比較から求められる．

$$\begin{array}{}\text{(9-24)}& \begin{array}{}\begin{array}{c}{k}_1=\underset{z\to \alpha }{\lim }\left(1-\alpha \mathrm{<mspace\; width="0.1em"></mspace>}{z}^{-1}\right)Y\left(z\right)={\left.\left(1-\alpha \mathrm{<mspace\; width="0.1em"></mspace>}{z}^{-1}\right)Y\left(z\right)\right|}_{z=\alpha }\\ k_2=\underset{z\to \beta }{\lim }\left(1-\beta {\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}\right)Y\left(z\right)={\left.\left(1-\beta {\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}\right)Y\left(z\right)\right|}_{z=\beta }\\ k_0=\frac{b_0}{a_0}-\left(k_1+k_2\right)\end{array}\end{array}\end{array}$$

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

$$\begin{array}{}\text{(9-25)}& y\mathrm{<mspace\; width="0.1em"></mspace>}\left[n\right]=k_0\delta \left[n\right]+k_1{\alpha }^n+k_2{\beta }^n\end{array}$$

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

$$\begin{array}{}\text{(9-26)}& Y\left(z\right)=\frac{k_1}{{\left(1-\alpha \mathrm{<mspace\; width="0.1em"></mspace>}{z}^{-1}\right)}^n}+\frac{k_2}{{\left(1-\alpha \mathrm{<mspace\; width="0.1em"></mspace>}{z}^{-1}\right)}^{n-1}}+\cdots +\frac{k_n}{1-\alpha \mathrm{<mspace\; width="0.1em"></mspace>}{z}^{-1}}\end{array}$$

ここで，k_r は以下の留数定理から求められる．

$$\begin{array}{}\text{(9-27)}& \begin{array}{}\begin{array}{r}{k}_r\\ \end{array}\begin{array}{l}=\frac{1}{\left(r-1\right)!}\underset{z\to \alpha }{\lim }\frac{d^{\left(r-1\right)}}{{dz}^{\left(r-1\right)}}{\left(1-\alpha \mathrm{<mspace\; width="0.1em"></mspace>}{z}^{-1}\right)}^{n}Y\left(z\right)\\ =\frac{1}{\left(r-1\right)!}\frac{d^{\left(r-1\right)}}{{dz}^{\left(r-1\right)}}{\left(1-\alpha \mathrm{<mspace\; width="0.1em"></mspace>}{z}^{-1}\right)}^n{\left.Y\left(z\right)\right|}_{z=\alpha }\end{array}\end{array}\end{array}$$

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

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

$$\begin{array}{}\text{(9-28)}& y\left[n\right]=\sum _{k=1}^N{a}_k\mathrm{<mspace\; width="0.1em"></mspace>}y\left[n-k\right]+\sum _{r=0}^M{b}_r\mathrm{<mspace\; width="0.1em"></mspace>}x\left[n-r\right]\end{array}$$

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

$$\begin{array}{}\text{(9-29)}& Y\left(z\right)=\sum _{k=1}^N{a}_k\mathrm{<mspace\; width="0.1em"></mspace>}Y\left(z\right)\mathrm{<mspace\; width="0.1em"></mspace>}{z}^{-k}+\sum _{r=0}^M{b}_r\mathrm{<mspace\; width="0.1em"></mspace>}X\left(z\right)\mathrm{<mspace\; width="0.1em"></mspace>}{z}^{-r}\end{array}$$

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

$$\begin{array}{}\text{(9-30)}& H\left(z\right)=\frac{Y\left(z\right)}{X\left(z\right)}=\frac{{\sum }_{r=0}^M{b}_r\mathrm{<mspace\; width="0.1em"></mspace>}{z}^{-r}}{1-\sum _{k=1}^N{a}_k{z}^{-k}}\end{array}$$

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

$$\begin{array}{}\text{(9-31)}& H\left(e^{j\omega T}\right)=\frac{{\sum }_{r=0}^M{b}_r\mathrm{<mspace\; width="0.1em"></mspace>}{e}^{-jr\omega T}}{1-\sum _{k=1}^N{a}_k{e}^{-jk\omega T}}\end{array}$$

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

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

$$\begin{array}{}\text{(9-32)}& G\left(s\right)=\frac{1}{\tau s+1}=\frac{{\omega }_c}{s+{\omega }_c}\end{array}$$

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

$$\begin{array}{}\text{(9-33)}& H\left(z\right)=\frac{{\omega }_c}{1-e^{-{\omega }_{c}T}{\mathrm{<mspace\; width="0.1em"></mspace>}z}^{-1}}\end{array}$$

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

$$\begin{array}{}\text{(9-34)}& H\left(e^{j\omega T}\right)=\frac{{\omega }_c}{1-e^{-{\omega }_{c}T}{e}^{-j\omega T}}\end{array}$$

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

図9-2: ディジタル化されたローパスフィルタのゲイン線図の例．

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

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

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

本章の参考文献

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