連続スペクトルの導出とその意味

連続スペクトルの導出とその意味
Z 変換の計算方法が分かったので、次は Z 変換を使って f [i] の連続スペクトル F(w) を求めてみます。
その後、連続スペクトル F(w) とはいったい何であるのか、その具体的な意味について考えます。
1
連続スペクトル F(w) の導出方法
非周期的で長さが無限であるデジタル時間領域信号 f [i] から Z 変換 F(z) を求めたら、それを元にして連続スペクト
ル F(w) を求めることができます。
連続スペクトルとは一体何なのかについては以前ラプラス変換の所で学びましたが、アナログ信号処理を履修してな
い人もいると思うので、導出の仕方を学んだら改めて詳しく学習することにします。
さて Z 変換の定義の所で話した様に「Z 変換はデジタル時間領域信号 f [i] を時間領域信号とみなしてラプラス変換し
たもの」なので、実は次のようにして簡単に Z 変換をラプラス変換に変形出来ます。
Z 変換 F(z) をラプラス変換 F(s) に変形
サンプリング間隔が τ (秒) の時、F(z) に z = esτ を代入して複素角周波数 s の関数にする。つまり
F(esτ )
を求める。これは f [i] のラプラス変換となっていて、改めて F(s) で表す。
なぜ z = esτ を代入するかについては説明を端折った所ですので、興味があれば自分で調べてみて下さい。
ともかく、これで f [i] のラプラス変換 F(s) が求まりましたので、後はラプラス変換の時と同様に s を jw に置き換
えればフーリエ変換、すなわち連続スペクトル F(w) が求まります。ここまでの話をまとめると
Z 変換 F(z) を連続スペクトル F(w) に変形
サンプリング間隔が τ (秒) の時、F(z) に z = ejwτ を代入して角周波数 w (rad/秒) の関数にする。つまり
F(ejwτ )
を求める。これは f [i] のフーリエ変換、すなわち連続スペクトルとなっていて、改めて F(w) で表す。
なお f [i] はデジタル信号であるのでナイキスト角周波数 ws /2 (rad/秒) より大きい角周波数のデジタルサイン波を含
まない。従って DFT の時と同様に F(w) も周期が ws で ws /2 を中心に複素共役関係 ( F(ws /2+x) = F∗ (ws /2−x), (0 ≤
x ≤ ws /2)) となる関数になる。よって w (rad/秒) の範囲は 0 ≤ w ≤ ws に限定して考える。
このように、「f [i] の w (rad/秒) における連続スペクトル値 F(w) は、 Z 平面上の点 z = ejwτ を F(z) に代入し
て求める」訳ですので、もし z = ejwτ が F(z) の収束領域に含まれていないと F(w) は発散します。この z = ejwτ の
位置についてはこの後の例題や演習の所で詳しく学びます。
逆に 0 ≤ w ≤ ws の範囲に含まれる全ての w (rad/秒) に関して連続スペクトル F(w) を求める事が出来るなら、も
のすごく大きい (でも有限な) ある正整数 N を i が超えると f [i] はほぼ 0 であるとみなす事が出来て (補足 1) 、かつ
F(w) は連続関数になります (補足 2) 。
1
2
DFT 係数 C[k] と 連続スペクトル F(w) の関係
連続スペクトル F(w) の意味を考える前に、DFT 係数 C[k] と 連続スペクトル F(w) の関係を見てみましょう。
なお f [i] は 0 ≤ w ≤ ws の範囲に含まれる全ての w (rad/秒) に関して連続スペクトル F(w) を求められることにし
ます。
ではまず C[k] 側から話を始めましょう。十分大きい正の整数 N を決めて f [0] から f [N − 1] まで切り出します。その
後、切り出した信号は周期 N の周期性デジタル信号であると仮定して N 点 DFT を行うと DFT 係数
C[k] =
N−1
∑
{f [i] · e−j
2π
N ·ki
}, (k = 0, 1, · · · , N − 1)
i=0
が求まります。
C[k] の話はここで一旦置いといて、次は F(w) 側の話に移ります。まずサンプリング角周波数 ws を N 等分し、角周
波数軸上に N 個の座標
wk =
ws
· k, (k = 0, 1, · · · , N − 1)
N
を定義します。各座標 wk における連続スペクトル F(wk ) の値は
F(wk ) =
∞
∑
{f [i] · e−jwk τ i }
i=0
なのですが、前節の最後で話したように、N が十分大きいければ f [i], (i ≥N) はほぼ 0 となりますので、思い切って
f [i] = 0, (i ≥N) に置き換えてしまいます。すると
F(wk ) =
N−1
∑
{f [i] · e−jwk τ i }
i=0
となります。ところでサンプリング間隔 τ (秒) は τ = 1/fs = 2π/ws より
wk =
2π
ws
·k =
·k
N
Nτ
と変形できますから、これを F(wk ) に代入して
F(wk ) =
N−1
∑
{f [i] · e−j
2π
N ·ki
}
i=0
となります。ここまで来たらもう分かりますね。つまり N が十分大きいければ
C[k] = F(wk )
です。ここまでの話をまとめましょう。
2
DFT 係数 C[k] と連続スペクトル F(w) の関係
f [i] が 0 ≤ w ≤ ws の範囲に含まれる全ての w (rad/秒) に関して連続スペクトル F(w) を求める事が出来るとす
る。また十分大きい正の整数 N を決める。
次に f [0] から f [N − 1] まで信号を切り出して N 点 DFT を行い、N 個の DFT 係数 C[k], (k = 0, 1, · · · , N − 1)
を求める。
一方、サンプリング角周波数 ws を N 等分し、角周波数軸上に N 個の座標
wk =
ws
· k, (k = 0, 1, · · · , N − 1)
N
を定義すると、
C[k] = F(wk ) =
N−1
∑
{f [i] · e−j
2π
N ·ki
}
i=0
が成り立つ。
この関係を例示したのが図 1 です (分かりやすいように 0 ≤ w ≤ ws /2 かつ連続振幅スペクトルだけ示しています)。
これを見ると分かるように
C[k] は F(w) をサンプリング間隔 ws /N (rad/秒) でサンプリングしたデジタル信号
となっています。
|F(w)|
7
|C[0]| = |F(0)|
6
|C[1]|= |F(w1)|
5
|C[2]| = |F(w2)|
4
3
|C[3]| = |F(w3)|
|C[5]| = |F(w5)|
2
1
|C[4]| = |F(w4)|
0
0
5
15
ws/N
ws/N
w0
10
w1
20
ws/N
ws/N
w3
w2
25
30
w
ws/N
w4
w5
図 1: DFT 係数 C[k] と 連続スペクトル F(w) の関係 (※)
(※) f [i] = {1, 2, 3, 0, 0, 0, · · · }、N = 10、fs = 10、τ = 0.1、 ws = 20π 、ナイキスト周波数 ws /2 = 10π 、ws /N = 2π
の条件で 0 ≤ w ≤ ws /2 における連続振幅スペクトル |F(w)| のグラフと |C[k]| の値を表示
3
連続スペクトル F(w) の意味
3
前節で話したように、十分大きい長さ N の DFT によって求めた C[k] は F(w) をサンプリング間隔 ws /N (rad/秒)
でサンプリングしたデジタル信号であり、 C[k] = F(wk ) の関係が成り立ちます
さてここで N を更に大きくして ∞ に近づけて行きます。すると座標 wk の間隔はどんどん狭まって 0 に近づいて行
きますので、うまく N と k を選んであげれば 0 ≤ w ≤ ws の範囲に含まれるどんな角周波数 w (rad/秒) の値でも wk
を使って表すことが出来るようになります。よって任意の角周波数 w (rad/秒) の F(w) を DFT 係数 C[k]) を使って表
せるようになります (補足 3) 。
ところで DFT 係数 C[k] は対応する角周波数のデジタルサイン波の振幅と位相を表す数字でしたね。以上から次のこ
とが言えます。
連続スペクトル F(w) の意味
f [i] が 0 ≤ w ≤ ws の範囲に含まれる全ての w (rad/秒) に関して連続スペクトル F(w) を求める事が出来るとす
る。また十分大きい正の整数 N を決める。
この時、f [i] には直流成分 F(0) と 0 < w ≤ ws /2(ナイキスト周波数) の範囲の全ての角周波数 w (rad/秒) のデジ
タルサイン波が含まれている。
具体的には角周波数 w (rad/秒) のデジタルサイン波は、連続振幅スペクトル |F(w)| と連続位相スペクトル ∠F(w)
を使って
gw [i] =
(
2|F(w)|
π)
sin w · τ · i + ∠F(w) +
, (i = 0, 1, · · · )
N
2
で表される。つまり「連続スペクトル F(w) は f [i] に含まれている角周波数 w (rad/秒) のデジタルサイン波の
振幅と位相を表している」。
ところで上で書いたデジタルサイン波 gw [i] の振幅の分母に N が入っているのに気が付きましたか? N は無限大に近
い大きな整数ですので、ひとつひとつの gw [i] の振幅はほとんど 0 になります。ただし w は連続量なので gw [i] は無限
に存在します。つまり振幅がほとんど 0 のデジタルサイン波が無限に足し合わされて f [i] が出来ているというイメージ
になります。
連続スペクトルの例
4
4.1
有限長 L の f [i] の連続スペクトル
長さが有限の L、つまり f [0] から f [L − 1] まで値が入っている (残りは f [i] = 0, (i = L, L + 1, · · · )) デジタル時間領
域信号の連続スペクトルは以前求めた F(z) にそのまま z = ejwτ を代入して
F(w) =
L−1
∑
{f [i] · e−jwτ i }
i=0
で求まります。
話が前後しますが、z = ejwτ が収束領域に含まれてない時は F(w) は求められませんでした。しかし有限長 L の f [i]
の場合は収束領域で無いのは原点 z = 0 だけでしたので、必ず z = ejwτ は収束領域に含まれている事が分かります ( x
が実数なら |ejx | = 1 でしたね)。
さて連続振幅スペクトルと連続位相スペクトルを求めるために z = e−jwτ i をオイラー公式を使って実数成分と虚数成
分に分けます。
4
F(w) =
L−1
∑
{f [i] · (cos(wτ i) − j sin(wτ i))} =
L−1
∑
[
{f [i] · cos(wτ i)} + j −
i=0
i=0
L−1
∑
]
{f [i] · sin(wτ i))}
i=0
これで実数成分と虚数成分が求まりましたので、後はいつも通り自乗和のルートで振幅スペクトル、atan2 関数で位相ス
ペクトルを求められます。
4.2
デジタルインパルス信号 f [i] = δ[i − n] の連続スペクトル
ある 0 以上の整数 n に対し、デジタルインパルス信号

1 (i = n)
f [i] = δ[i − n] =
0 (i 6= n)
の連続スペクトルは求めてみましょう。n によって場合分けします。
・ n = 0 の場合
F(w) = 1(= |1| · e0 )
なので、振幅スペクトルは常に |F(w)| = 1、位相スペクトルは常に ∠F(w) = 0 です。
・ n 6= 0 の場合
F(w) = e−jwτ n = |1| · ej(−wτ n) )
なので、振幅スペクトルは常に |F(w)| = 1、位相スペクトルは ∠F(w) = −wτ n です。
4.3
等比数列 f [i] = ai の連続スペクトル
等比数列 f [i] = ai の Z 変換の収束領域は前 2 つの例と違って少し面倒な形になっていましたので、連続スペクトルの
収束についても真面目に考える必要があります。
等比数列 f [i] = ai の Z 変換の収束領域は前のテキストに示したように Z 平面上の半径 |a| の円の外側でしたので、
z = ejwτ も半径 |a| の円の外側にある必要があります。ところで w を 0 ≤ w ≤ ws の範囲で変化させると z = ejwτ は
z = 1 からスタートして半径 1 の単位円の上を半時計回りに回転する円運動になります。従って「等比数列 f [i] = ai は
|a| < 1 ならば連続スペクトルが存在して、 |a| ≥ 1 ならば存在しない」ことが分かります。もっとも、わざわざ収束領
域とか考えなくても、 |a| = 1 なら f [i] は無限に 1 が出続ける数列になりますし、|a| > 1 なら f [i] は無限大に発散す
るので、連続スペクトルが無いのは当然ではありますが。
さて |a| < 1 のとき連続スペクトルが存在して、その値は収束後の F(z) の式に z = ejwτ を代入して
F(w) =
1
1 − a · e−jwτ
となります。この後は最初の例と同様に連続振幅スペクトルと連続位相スペクトルを求めるために実数成分と虚数成分
に分けるのですが、今回の Z 変換は分母に z 以外の項が入っているので少し計算が面倒です (補足 4) 。結論を書くと
F(w) =
1 − a cos(wτ )
−a sin(wτ )
+j·
1 − 2a cos(wτ ) + a2
1 − 2a cos(wτ ) + a2
です。
5
演習
演習 1 (個人):C 言語を使って、有限長 L = 3、f [i] = {1, 2, 3, 0, 0, · · · } の連続振幅スペクトル |F(w)| を戻り値とする
関数 double abs F( double w ) と連続位相スペクトル ∠F(w) を戻り値とする関数 double arg F( double w ) を作成せ
よ。なお L と τ = 0.1 は外部変数とする。
更にそれらの値を printf で出力し (フォーマットは”w, 振幅, 位相¥n”。csv ファイルに出力すればよりベター)、連続
振幅スペクトルと連続位相スペクトルのグラフを表計算ソフトを使って描いて保存せよ。w の範囲は 0 ≤ w ≤ ws /2 と
し、代表点は 1 (rad/秒) 刻みとする (まずはナイキスト角周波数 ws /2 がいくつになるかから考えよう)。
演習 2 (チーム):(前問の続き) DCT 係数が連続スペクトルのデジタル信号となっていることを視覚的に確認する。L = 3、
f [i] = {1, 2, 3, 0, 0, · · · } の場合は DFT 点数 N としていくつ以上に設定すれば十分であるかチームで考えてノートにま
とめよ。
演習 3 (個人):(前問の続き) C 言語を使って L = 3、f [i] = {1, 2, 3, 0, 0, · · · } の N 点 DFT を求め、k = 0, 1, · · · , N/2
に対して振幅成分と位相成分を求めよ。同時に wk =
ws
N
· k も求めよ。
更にそれらの値を printf で出力し (フォーマットは”wk , 振幅, 位相¥n”。csv ファイルに出力すればよりベター)、振
幅成分と位相成分のグラフを表計算ソフトを使って描いて保存せよ。
演習 4 (個人):表計算ソフトを使って a = 0.1、f [i] = (0.1)i のグラフを描いて保存せよ。i の範囲は i = 0, 1, · · · , 100
とする。
演習 5 (個人):(前問の続き) a = 0.1、f [i] = (0.1)i の収束領域をノートに書け。また w = 0, ws /8, ws /4, 3ws /8, ws /2, 3ws /4, ws
の時の z = ejwτ の位置を点で表わし、それらの点が収束領域内にある (よって Z 変換が収束する) ことを確認せよ。
演習 6 (個人):(前問の続き) C 言語を使って、a = 0.1、f [i] = (0.1)i , (i = 0, 1, · · · ) の連続振幅スペクトル |F(w)| を戻
り値とする関数 double abs F( double w ) と連続位相スペクトル ∠F(w) を戻り値とする関数 double arg F( double w )
を作成せよ。なお a と τ = 0.1 は外部変数とする。
更にそれらの値を printf で出力し (フォーマットは”w, 振幅, 位相¥n”。csvァイルに出力すればよりベター)、連続振
幅スペクトルと連続位相スペクトルのグラフを表計算ソフトを使って描いて保存せよ。w の範囲は 0 ≤ w ≤ ws /2 とし、
代表点は 0.01 (rad/秒) 刻みとする。
演習 7 (チーム):(前問の続き) DCT 係数が連続スペクトルのデジタル信号となっていることを視覚的に確認する。
a = 0.1、f [i] = (0.1)i の場合は DFT 点数 N としていくつ以上に設定すれば十分であるか「3 つ前の演習で描いた f [i]
のグラフを見ながら」チームで考えてノートにまとめよ。
演習 8 (個人):(前問の続き) C 言語を使って a = 0.1、f [i] = (0.1)i の N 点 DFT を求め、k = 0, 1, · · · , N/2 に対して
振幅成分と位相成分を求めよ。同時に wk =
ws
N
· k も求めよ。
更にそれらの値を printf で出力し (フォーマットは”wk , 振幅, 位相¥n”。csv ファイルに出力すればよりベター)、振
幅成分と位相成分のグラフを表計算ソフトを使って描いて保存せよ。
演習 9 (個人):(前問の続き) 「N の値は前問のまま」としながら、a を a = 0.99 とした時の f [i] 、連続スペクトル、
DCT 係数のグラフを表計算ソフトを使って描いて保存せよ。
演習 10 (チーム):以上の演習で作成した 3 パターンの連続スペクトルと DCT 係数のグラフを見比べて、何が言えるの
かチームで話し合ってノートにまとめよ。
6
補足
本題から外れるので本文中では詳しくは説明しなかった事項を補足します。興味があったり、大学編入試験対策でもっ
と深く知りたい学生は読んで下さい。
補足 1:f [i] が 0 ≤ w ≤ ws の範囲に含まれる全ての w (rad/秒) に関して連続スペクトル F(w) を求める事が出来ると
いうのは、f [i] が次の「絶対収束条件」を満たしていることを意味します。
∞
∑
|f [i]| < ∞
i=0
それで、「f [i] が絶対収束するなら級数も収束する」という定理を使うと
∞
∑
f [i] < ∞
i=0
次に、「f [i] の級数が収束するなら f [i] は 0 に収束する」という定理を使うと
f [i] → 0, (i → ∞)
このことは -δ 論法で言い換えると、「もう 0 と見なしても良いよねという小さな数字 > 0 をユーザーが事前に決め
た時、|f [i]| < , (i ≥ N) となる整数 N が必ずある」ことを意味します。この N がいくつになるかは をある人がどう
選ぶかによって変わりますので、同じ f [i] でも人によっては N = 100 程度で済む場合もあるし、人によっては N = 全
宇宙の星の数とかになったりします。
補足 2:f [i] が 0 ≤ w ≤ ws の範囲に含まれる全ての w (rad/秒) に関して連続スペクトル F(w) を求める事が出来ると
き、f [i] は絶対収束条件
∞
∑
|f [i]| < ∞
i=0
を満たしていると一つ前の補足で言いましたが、これをフーリエ変換の用語で言い換えると「f [i] をアナログ信号と
みなした時 f [i] は絶対可積分である」という事を意味します。
一方「ある関数が絶対可積分であるならそのフーリエ変換は一様連続だ」という定理があるので、F(w) は連続関数だ
という事が分かります。
補足 3:一つ前の補足で f [i] が 0 ≤ w ≤ ws の範囲に含まれる全ての w (rad/秒) に関して連続スペクトル F(w) を求め
る事が出来るなら F(w) は連続関数になると言いました。
一方 N を無限大に近づけると wk 全体の集合は 0 ≤ w ≤ ws の稠密な部分集合になります。稠密というのは 0 ≤ w ≤ ws
の範囲のどんな w でも wk をつかって無限の精度で近似できるという意味です。
ここで「ある連続関数の任意の値は定義域の稠密な部分集合上の関数値のみで決まる」という定理を使うと、どんな
w でも F(w) の値は F(wk ) 、つまり C[k] を使って表すことが出来ます。
補足 4:次のように F(w) = 1/(1 − ae−jwτ ) の分母と分子に 1 − aejwτ をかけて分母から虚数成分を取り除いてからオ
イラー公式を適用します。
1
1 − aejwτ
1
=
·
−jwτ
−jwτ
1 − ae
1 − ae
1 − aejwτ
jwτ
1 − a cos(wτ ) − j · a sin(wτ )
1 − ae
=
=
jwτ
−jwτ
2
1 − a(e
+e
)+a
1 − 2a cos(wτ ) + a2
F(w) =
7