信号処理ノート
静岡大学工学研究科
機械工学専攻
大坪 順次
1.繰り返しの信号とフーリエ級数展開 1.1 繰り返しのある信号 1.2 信号の複素表示 1.3 複素表示を用いた計算の例 1.4 フーリエ級数展開 1.5 指数関数を用いた級数展開 1.6 フーリエ級数展開の例 2.デルタ関数とフーリエ変換 2.1 デルタ関数 2.2 デルタ関数を使った関数値の計算 2.3 デルタ関数のさまざまな定義 2.4 フーリエ変換 2.5 フーリエ変換の例 2.6 フーリエ変換の性質 2.7 離散フーリエ変換 2.8 高速フーリエ変換−FFT 2.9 FFT の実際の計算方法 3.線形システムと伝達関数 3.1 線形システム 3.2 線形システムにおけるたたみ込み積分 3.3 たたみ込み積分の意味 3.4 伝達関数 4.信号の相関、パワースペクトル、サンプリングの定理 4.1 信号相関 4.2 相関関数の例 4.3 相関関数とパワースペクトル 4.4 信号のサンプリング 4.5 サンプル関数のフーリエ変換
4.6 サンプリングの定理 4.7 サンプリングの定理の意味 5.信号の再生、回復 5.1 デコンボリューション 5.2 反復法を用いたデコンボリューション 5.3 逆フィルタ 5.4 ウィーナーフィルタ 5.5 修正フィルタ 6.不規則信号 6.1 雑音と確率過程 6.2 エルゴード性 6.3 平均値と分散 6.4 二変数相関と共分散 6.5 正規分布 6.6 雑音の重ね合わせ 6.7 実際の雑音とSN 比 7.雑音除去 7.1 雑音除去の方法 7.2 雑音の積算平均化 7.3 雑音除去フィルタ 7.4 周波数領域での雑音除去 7.5 信号波形の検出
1 繰り返しの信号とフーリエ級数展開 1.1 繰り返しのある信号 計測やデータ処理において、さまざまな信号が発生し、また加工される。こ れらの信号の振幅はsin あるいは cos 関数などの重ね合わせとして表すことが可 能である。正弦関数、余弦関数は波の伝搬などにおいても頻繁に使われている。 ここでは、信号を一つの波として考え、まず波の伝搬について一般的な性質と その表記から入ろう。 ある座標x 軸方向にv なる速度で進む波は € ψ(x,t) = f (x − vt) (1.1) のように書ける。t=0 は、空間的波の波形を表す。(1.1)式を満たす波の例として、 次のような余弦波が考えられる。 € ψ(x,t) = A0cos{k(x − vt) +φ} = A0cos(kx −ωt +φ) (1.2) ここで、A0は振幅の大きさ、v=ω/k、φは波の初期位相、ωは角周波数、k は波数 を表す。1.3 節で示すように、繰り返しのある任意の波形の波として、余弦関数 あるいは正弦関数で表される波の重ね合せとして表すことができる。すなわち、 (1.2)式のような余弦波と正弦波の線形和として一般の波が表されるということ である。波は、余弦、正弦関数以外の直交関数のセットを用いても直交展開す ることができる。ここで正弦、余弦関数を使うのは、波を合成するのにもっと よく使われる関数だからである。波は直交関数の線形和で表されるから、その 成分の一つである(1.2)式の波についていろいろなことを調べると、一般の波に ついて調べたことと同じであると言える。したがって、当面は波といったとき に(1.2)式で表される波の性質について調べればよいことにする。 1.2 信号の複素表示 ここで、波を表すのに(1.1)式の余弦関数を用いたが、余弦関数の代りに複素 数で波を表示し、その実部が物理的な波動を表すものとすると、いろいろな波 動の振舞を楽に計算することができる。ここでは、しばらくの間波動の複素数 表示について考えてみよう。図1.1 に示すように直交座標において、あるベクト ルc は複素表示を用いて €
c = a + ib, a = Re[c], b = Im[c] (1.3)
€ r = a2 + b2 (1.4a) € θ = tan−1b a (1.4b) 図1.1 フェザー表示 のように表せる。ここで、cosθ=a/r、sinθ=b/r である。また、r は c の絶対値 € | c |= r、 θは偏角 arg(c)=θである。オイラー(Euler)の公式を用いると € c = r exp(iθ) (1.5a) €
sinθ =exp(iθ) − exp(−iθ) 2i = Im c r ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ (1.5b) €
cosθ =exp(iθ) + exp(−iθ) 2 = Re c r ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ (1.5c) となる。(1.5a)式のような表示法をフエザー表示という。(1.5a)式の指数部を負の 符号で置換えたものをc の複素共役と言い、c*で表す。すなはち €
c* = r exp(−iθ) = r{cos(−θ) + isin(−θ)} = a − ib (1.6) である。フエザー表示によると、(1.2)式は、A0を実数として € ψ = Re[A0exp{i(kx −ωt +φ)}] (1.7) と書ける。A= A0exp{i(kx-ωt+φ)}を波動の複素振幅と言う。 波動の足し算などは複素数として計算することができるが、計算ができるの はあくまでも線形演算に限られる。たとえば、波動の平均パワーを求めるとき にはどのようになるだろうか。実関数の場合には、波が(1.2)式のような cos 関数 で表されるものとして、その時間平均パワーは、波の周期に渡っての積分平均、 すなわち
€ I =< ψ2 >= ω 2π A0 2 cos2 −π /ω π /ω
∫
(kx − ωt + φ)dt = A0 2 (1.8) で計算される。一方、複素数で計算する場合には、単に複素表示された波動を2 乗して平均値を計算したのでは、(1.11)式と同じ結果にならない。そのため、(1.2) 式の波をA として複素表示したとき、便宜上そのパワーを € I = AA * 2 = A02 2 (1.9) と対応させることにする。 1.3 複素表示を用いた計算の例 波の重ね合せについて考えてみよう。一般的には、異なる波形の波も(1.5a)式 に示すようなフエザーの加算として表される。ここでは、簡単のため同じ波形 ではあるが波の大きさ A0 と位相φが異なる次の二つの波同士の重ね合せについ て考えてみよう。二つの波 €ψ1 = A01cos(kx − ωt + φ1) = Re[A01exp{i(kx − ωt) + iφ1}] (1.10a)
€
ψ2 = A02cos(kx − ωt + φ2) = Re[A02exp{i(kx − ωt) + iφ2}] (1.10b) の重ね合わせは、複素振幅を用いて計算すると
€
ψ =ψ1+ψ2 = Re[A1+ A2]
=| A01exp(iφ1) + A02exp(iφ2) | cos{kx −ωt + arg(A1) − arg(A2)}
(1.11) となる。余弦関数のまま波動の重ね合わせを計算することは手間がかかること が多いが、フエザー表示を使うと計算が簡単にできる利点がある。 実 際 に 、 便 利 に 計 算 で き る 具 体 的 な 例 題 を 解 い て み よ う 。 二 つ の 波 € ψ1 = 4 cos(kx −ωt)と € ψ2 = 3cos(kx −ωt +π/2)となる波の和を考える。cos 関数のま ま波を合成しても計算はできるが、三角関数の公式を思い出しながらの計算に なり、また計算間違いもしやすい。一方、これを複素振幅として解くと €
ψ =ψ1+ψ2 = Re[4 exp{i(kx −ωt)} + 3exp{i(kx −ωt) + i π
2}] = Re[| 4 + i3 | exp{i(kx −ωt + tan−13
4)}] = 5cos(kx −ωt + tan −13 4) (1.12) と、非常に簡単に答えを導くことができる。この足し算の関係を、複素平面で 表すと、図1.2 のようになる。信号処理においては、信号、関数を複素表示とし て取り扱うことが多い。複素表示を使うことの利点としては、線形操作におい
て実関数(sin、cos などの三角関数)の和や差、あるいは積分微分を計算が多く 出てくる。その場合に、複素表示を用いて計算することによって計算を簡略化 できるとともに、計算量も激減させることができる。このとき、実際の物理的 信号としては、計算した関数の実部を考えればよい。 図1.2 複素計算の例 1.4 フーリエ級数展開 以下では、周期的関数のフーリエ級数展開について述べる。ある関数f(x)が間 隔x0でくり返す周期性の波であるとき、f(x+x0)=f(x)のように書くことができる。 このような関数f(x)は周期関数と呼ばれる。このとき、次のような関数は一般に 三角関数の無限和で表すことができる。 € f (x) = a0 x0 + 2 x0 n =1{ancos(n2πνx) + bnsin(n2πνx)} ∞
∑
(1.13) ここで、ν=1/x0である。この展開は、直交展開の正規直交関数のセットとして、 cos、sin 関数を選んだことに他ならない。このような展開をフーリエ級数展開と いう。三角関数の直交性を使うと、展開係数an、bnは € an = −x f (x)cos(n2πνx)dx 0/ 2 x0/ 2∫
(n=0,1,2,…) (1.14a) € bn = −x f (x)sin(n2πνx)dx 0/ 2 x0/ 2∫
(n=1,2,3,…) (1.14b) として求めることができる。一般的に、系の物理的な応答を調べるときに、し ばしば単一のみの関数cos あるいは sin 関数で表される振動子の応答についての み記述されることが多い。このことは、(1.13)式からわかるように、任意の周期 関数は、基本周波数をνとして、その高周波成分との調和振動子の線形結合で表されるという事実に基づいている。このことにより、ある特定の周波数の調和 振動子について、その系における振舞いを調べておけば、それらの線形結合と してその物理現象が理解されるというものである。 1.5 指数関数を用いた級数展開 (1.13)式の cos、sin 関数の代わりに、オイラーの関係式を使って、指数関数を 用いてフーリエ展開を書き表してみよう。cos、sin 関数を指数関数で表して、 (1.13)式に代入すると、f(x)は € f (x) = a0 x0 + 1 x0 (an − ibn)exp(i 2πnx x0 n =1 ∞
∑
) + 1 x0 (an+ ibn)exp(−i 2πnx x0 n =1 ∞∑
) (1.15) のように変形できる。ここで A0 = a0 (1.16a) € An = an − ibn (1.16b) € A−n = an + ibn (1.16c) のように、複素フーリエ展開における係数を定義すると € f (x) = 1 x0 Anexp(i 2πnx x0 n =−∞ ∞∑
) (1.17) と書くことができる。一般に、くり返しの関数は、n が−∞から∞にわたる複素 係数 An(これをスペクトルという)を持つ指数関数からなる振動項の重ね合わ せによって展開できることを示している。このことは、指数関数も正規直交関 数のセットとなっていることを示している。関数f(x)が与えられているとき、係 数Anは € An = f (x)exp −i2nπx x0 ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ dx −x0/ 2 x0/ 2∫
(1.18) で与えられる。(1.17)式は複素数での展開となっているが、A-n=An*の関係が成り 立つときには、f(x)は実関数となることがわかる。この式は、式(1.13)と等価であ るが、実関数の波動を複素振幅として表したときと同様、(1.17)式を用いた方が 計算が楽になる利点があり、通常の繰り返しのある波の計算においては、(1.17) がよく使われる。 1.6 フーリエ級数展開の例 フーリエ級数展開の例として、周期x0、デューティサイクル1/2(パルス間隔のハイとローの比が同じ)である矩形波をフーリエ級数展開で表してみよう。 図1.3 矩形波 図1.3 に示す矩形波の定義を € f (x) = 1 nx0 ≤ x < (n +1 2)x0 −1 (n +1 2)x0≤ x < (n +1)x0 ⎧ ⎨ ⎪ ⎩ ⎪ (1.19) とする。(1.17)、(1.18)式の定義にしたがって、この関数のフーリエ展開は € f (x) = 2 iπ 1 2n −1exp{i2π(2n −1) x x0} n =−∞ ∞
∑
(1.20) と表される。(1.19)式の定義では、区分としてしか関数が定義できないが、(1.20) 式では、級数展開の形ではあるが、一つの関数として区分をすることなく定義 ができている。図1.4(a)は、(1.20)式を実数式(1.13)の形式で € f (x) = 4 π 1 2n −1sin{2π(2n −1) x x0} n =1 ∞∑
(1.21) と表し、級数和の各成分の波形をn について 1 から 4 までを表したものである。 また、図1.4(b)は、n について 1 から 5 までの成分を加算した結果を表している。 n=5 までの和ではまだ十分ではないが、矩形波の形が見える。加算する級数和の 項を増やすにつれ、合成波形は図1.3 に示した矩形波に漸近していく。 図1.4 矩形波の分解(a)と合成(b) f(x) x x0/2 x0 0 1 −1(1.13)式あるいは(1.17)式を見ると、くり返しの波動は、くり返しの基本周波数 1/x0 とその整数倍の高次の波を重ね合わせにより表されていることがわかる。 (1.17 )式の展開係数 an、bnや(1.17)式の Anは関数f(x)のスペクトルと呼ばれ、調 和振動成分がどの程度その関数に含まれているかを示す指標になっている。図 1.5 は、(1.20)式で表される関数のスペクトルを示したものである。くり返しの 波は、このように図1.5 に示すようなスペクトルを重みとして、基本波とその高 次の調和振動の和として表される。あるいは、くり返しの波は、基本波とその 高次振動に分解して考えることができるという言い方もできる。 図1.5 矩形波のスペクトル 1.7 関数の直交展開 波は、たとえば cos、sin 関数の合成として表すことができることを示した。 話は前後するが、波に限らず任意の関数に関して、このような直交関数による 展開ができるということについて、簡単に説明しておこう。あるなめらかな関 数f(x)は、関数が定義される空間における正規直交関数のセット{ψn}を使って € f (x) = Fnψn(x) n =−∞ ∞
∑
(1.22) のように直交展開することができる。Fn は展開係数である。また、直交関数の 定義より € ψn(x)ψn * (x)dx =δnm −∞ ∞∫
(1.23) である。ここで、δnmはクロネッカー(Kronecker)のデルタであり、n=m のときδnm=1 で、 ではδnm=0 である。f(x)に同じ次数の直交複素共役関数をかけ積分する ことによって、Fnは€ Fn = f (ξ)ψn * (ξ)dξ −∞ ∞
∫
(1.24) として求められる。この関係は物理的な意味合いも持っているが、ここでは一 般的な数学の関係として簡単にふれるだけにしておこう。 ところで、先に述べた繰り返しの関数において、sin、cos を用いて展開したと いうことは、繰り返しの関数が定義される空間で、sin、cos 関数が直交系となっ ているという意味である。(1.14)式の関係は、まさに(1.23)式と同じである。同様 に、複素指数関数においても exp 関数が直交系になっているということを意味 している。(1.17)式は、(1.13)式を異なる直交表現で表していることに他ならない。 したがって、一元の繰り返しの波の定義される空間においては、これ以外にも さまざまな直交系を導入することが出来る。また、新しく導入された直交系は、 数学の定理により他の直交系との線形変換が可能である。sin、cos 関数、あるい は複素指数関数を選択したのは、単に後々の計算がやりやすいようにと導入し たに過ぎない。他の目的で計算が容易になるようにするためには、新しい直交 関数を改めて導入することも出来る。 もう一つ、よく教科書において用いられる記述について述べておこう。我々 が目にする波、信号は正弦波関数のような単調な信号ばかりではなく、繰り返 しの波であっても先に見た矩形波のように単純なsin、cos の信号とは限らない。 しかし、波の伝搬、信号の伝搬の取り扱いにおいては、一般的な波、信号の形、 関数をそのまま仮定せずに、sin 関数あるいは cos 関数として信号を仮定した議 論が用いられる。一般の信号はそのような簡単な波形ではないのに、なぜその ような議論をするかというのは以下のとおりである。(1.13)式からわかるように、 一般の繰り返しのある信号は、sin、cos の信号の線形重ね合わせ(直交展開)と して表示することができる。したがって、その要素の一つの直交関数について 議論すれば、その線形和である関数f(x)について議論したことと同じなるという のがその理由である。このことから、波の伝搬、信号処理においては、その信 号を構成するある特定の周波数についての性質を知るというのが非常に重要な ポイントになる。 課題1 図 1.3 の矩形波(1.19)式をフーリエ級数展開して、(1.21)式で表すことが できる。このことを計算により示せ。また、複素表示として、(1.20)式のように2.デルタ関数とフーリエ変換 2.1 デルタ関数 本節では、フーリエ級数展開をフーリエ変換に拡張するが、その前にデルタ 関数について述べる。フーリエ変換、フーリエ光学では、デルタ関数がさまざ まな計算過程において出てくる。このデルタ関数は、フーリエ変換にとどまら ず応用範囲が広く、便利な関数である。デルタ関数の定義としては、ある唯一 の特別な関数形が与えられているわけではなく、下記の関係を満たすすべての 形式がデルタ関数であると定義される。 € δ(x) = 0 (x ≠ 0) δ(x)dx −∞ ∞
∫
= 1 ⎫ ⎬ ⎪ ⎭ ⎪ (2.1) したがって、デルタ関数は解析的な関数というよりは、ある関数として x=0 の 近傍でのみで零でない値をとり、x=0 において特異点となる超幾何関数として定 義される。厳密にいうと、式(2.1)は概念的なデルタ関数の定義であり、この ような超関数の定義として十分とはいえず、以下に述べるように、積分の中で はじめて意味のある定義がされる関数である。このような関数の候補はたくさ んある。たとえば、制御工学において広く用いられているステップ関数を微分 したものを考えると € du(x) dx = 0 (x ≠ 0) (2.2a) € du(x) dx dx −∞ ∞∫
= u(x)[
]
−∞ ∞ = 1 (2.2b) であるから、 (2.1)式の定義に従えば、これはデルタ関数であるといえる。デル タ関数の最も簡単な定義としては、次の形式がしばしば使われる。すなわち、 微小な量ε( € ε →0)に対し € δ(x) = 1 2ε | x |≤ε 0 | x |>ε ⎧ ⎨ ⎪ ⎩ ⎪ (2.3) で定義される関数は、(2.1)式の定義よりデルタ関数である。この定義によると、 デルタ関数は図2.1 に示すような関数形である。2.2 デルタ関数を使った関数値の計算 次に、デルタ関数の最も有用な性質を述べておこう。f(x)を微分可能かつ連続 な解析的関数とすると、 に対して (2.3)式を使い € f (x)δ(x)dx −∞ ∞
∫
= 1 2ε −ε f (x)dx = f (θ) → f (0) (ε →0) ε∫
(2.4) となる関係が得られる。ここで、εは微小量として、区間[−ε,ε]で関数 f(x)を線形 近似できるものとし、f(θ)={f(ε)+f(−ε)}/2 と仮定した。 (2.4)式から、ある関数と デルタ関数とを掛け合わせて積分したものは、元の関数の x=0 のときの値に等 しいという結果が得られる。このことから、さらに次に述べる重要な結果が得 られる。すなわち、デルタ関数δ(x)の座標を a だけ移動させたものと関数 f(x)を かけ合わせ、 (2.4)式と同様な計算をすると € f (x)δ(x − a)dx = f (a) −∞ ∞∫
(2.5) となる。このように、座標シフトさせたデルタ関数と任意関数の (2.5)式の積分 を行うと、元の関数f(x)が再生されることがわかる。デルタ関数が偶関数である ことは、次のようにして容易に確かめられる。 図2.1 δ 関数 € f (x)δ(−x)dx = f (−x)δ(x)dx = f (0) −∞ ∞∫
−∞ ∞∫
(2.6) また、デルタ関数の表現として重要なものとして € δ(x) = exp(−i2πνx)dν −∞ ∞∫
(2.7) の形のものがある。これが、デルタ関数の定義になっていることは、つぎのよ うにして確かめられる。すでに述べたように、関数f(x)は、関数が定義される空間の正規直交関数のセット{ψn}を使って € f (x) = Fnψn(x) n =−∞ ∞
∑
(2.8) のように直交展開することができる。Fn は展開係数である。また、直交関数の 定義より € ψn(x)ψn * (x)dx =δnm −∞ ∞∫
(2.9) である。ここで、δnmはクロネッカー(Kronecker)のデルタであり、n=m のときδnm=1 であり、 € n ≠ mではδnm=0 である。Fnは、f(x)に同じ次数の直交複素共役関数をか け積分することによって € Fn = f (ξ)ψn * (ξ)dξ −∞ ∞∫
(2.10) として求まる。 これらのことを使うと、f(x)は形式的に € f (x) = f (ξ) ψn*(ξ)ψn(x)dξ n =−∞ ∞∑
−∞ ∞∫
(2.11) のように書くことができる。この式と (2.5)式を比べると、積分の中の級数和の 部分はデルタ関数であることがわかかる。すなわち € δ(ξ − x) = ψn *(ξ)ψ n(x) n =−∞ ∞∑
(2.12) である。ここで、無限区間で定義された直交関数として € ψn(x) = 1 2N exp(i nπ N x) (2.13) と選ぶと、デルタ関数は € δ(ξ − x) = 1 2N exp{−i nπ N (ξ − x)} n =−∞ ∞∑
(2.14) と与えられる。ここで、 と置くと、N を限りなく大きくすると、 (2.14) 式は € δ(ξ − x) = exp{−i2πνn(ξ − x)}Δνn n =−∞ ∞∑
→ exp{−i2πν(ξ − x)}dν −∞ ∞∫
(2.15) となり、確かに (2.7)式がデルタ関数の定義になっていることがわかかる。ここ で、 € Δνn =νn +1−νnとした。上記の証明は数学的な精密さを欠くが、 (2.7)式の結果は非常に有用である。 たとえば、後でも示すようにフーリエ変換などの計算において、この結果を使 い積分の次元を減らすことができる。また、この式の意味するところは、一定 の値の関数(定数)のフーリエ変換はデルタ関数となることを表している。実 時間信号が近似的にデルタ関数(実際にはt=0 で無限大に発散する信号なので、 そのような信号は実現できない)として与えられたとすると、そのフーリエ変 換である周波数スペクトル成分は、0 周波数から∞に至る周波数まで一様な振幅 で記述されることになる。実際には発散してしまう信号になってしまうが、高 い周波数成分(「高い」の意味はそれぞれ考えるシステムでの相対的なレベルに なるが)まで一様とみなされるスペクトルを持つ信号はデルタ関数に近似でき るといえる。また逆に、無限大時間にわたって一様なDC 信号(実際には無限時 間は実現できないし、長時間一定といってもドリフトが存在する)があれば、 その信号のスペクトルは 0 周波数だけの成分を持つデルタ関数で与えられるス ペクトルとなる。 2.3 デルタ関数のさまざまな定義 この他にも、デルタ関数として定義できる形がいろいろあるが、それらの多 くは、解析的な関数のある極限として定義される。以下に、1 次元のδ関数のい くつかの例をあげる。2 次元のδ関数についても、これを拡張して同様に定義す ることができる。 € δ(x) = N exp(−N2πx2) (N →∞) (2.16) € δ(x) = Nrect(Nx) = 1 | x |≤ 1 2N 0 | x |> 1 2N ⎧ ⎨ ⎪ ⎩ ⎪ (N →∞) (2.17) € δ(x) = Nsinc(Nx) = NsinπNx πNx (N →∞) (2.18)
ここで、rect(x)は rectangular 関数と呼ばれる。また、sinc(x)は sinc 関数と呼ばれ、
€ sinc(x) = sin(πx) /(πx)のように定義される関数である。これらの関数形を図 2.2 に示す。これらを含むデルタ関数近似は、積分計算などにおいて非常に重要か つ使い道がある。たとえば、多重積分において、上記のような関数が含まれる 場合に、これをデルタ関数と近似し、(2.5)式の関係を使うことにより積分の次 元を一つ少なくすることができ、積分計算をより簡単にすることができる。
図2.2 δ 関数の例 (a)ガウス関数近似、(b)矩形関数近似、(c)sinc 関数近似 2.4 フーリエ変換 次に、フーリエ級数展開を拡張し、関数のフーリエ変換について述べる。(1.13) 式の展開では有限な周期を持つ関数について考えたが、周期 x0が限りなく大き くなる。すなわち、関数が周期的ではなくなるとき、(1.14)式あるいは(1.18)式で 表されるスペクトルの隣会う間隔は、n/x0について限りなく小さくなる。すなわ ち、(1.18)式の Anはν=n/x0を(ω=2πν)振動数とする連続するスペクトル密度と考 えることができる。この連続するスペクトル密度をF(ν)で表すと、式(1.17)は € f (x) = F(ν)exp(i2πνx)dν −∞ ∞
∫
(2.19) のように表すことができる。デルタ関数の関係を使うと、f(x)が知られていると き、f(x)に exp(−i2πνx)をかけ[−∞,∞]で積分すると € f (x)exp(−i2πνx)dx = dx −∞ ∞∫
−∞ ∞∫
dν' F(ν')exp{i2π(ν' −ν)x} −∞ ∞∫
=∫
−∞∞F(ν')δ(ν' −ν)dν' = F(ν) (2.20) が得られ、スペクトル密度を計算することができる。(2.19)、(2.20)式のような関 係にある関数 f(x)と F(ν)とはフーリエ変換の関係にあるといい、(2.20)式を f(x) のフーリエ変換、(2.19)式を F(ν)の逆フーリエ変換と呼んでいる。フーリエ級数 はくり返しのある波を表すのに用いたが、フーリエ変換は無限大の周期の波、 すなわち、くり変えしのない単発現象の波(たとえば一つのインパルス波)の ような関数を表すのに用いることができる。 2.5 フーリエ変換の例 いま、図2.3(a)に示す[−a/2,a/2]の区間でのみ値を持つパルス波 € f (x) = 1 | x |≤ a /2 0 | x |> a /2 ⎧ ⎨ ⎩ (2.21)を考えてみよう。式(2.20)の関係を用いると € F(ν) = exp(−i2πνx)dx = asinc(aν) −a / 2 a / 2
∫
(2.22) なるスペクトル密度が得られる。F(ν)は、図 2.3(b)に示すような関数である。(2.21) 式のような単発波形の非周期的関数は連続スペクトルを持ち、スペクトルの振 幅はその振動数νにおける波の密度を表している。元々フーリエ変換は、ここで みたように単発の波のスペクトル構造を解析するために導入されたものである が、フーリエ級数展開できる周期的な関数であっても、これを適用し、そのス ペクトル解析を行うことも可能である。 2.6 フーリエ変換の性質 ここでは、フーリエ変換に関するいつくかの有用な性質について述べる。原 関数をf(x)とし、そのフーリエ変換を F(ν)とする。フーリエ変換された関数 F(ν) は、一般的に複素関数となり € F(ν) =| F(ν) | exp{iφ(ν)} (2.23) のように書ける。原関数が実数で表される物理的な内容を持つ実関数であって も、そのフーリエ変換は一般に複素関数である。したがって、φ(ν)は複素関数 F(ν)の位相成分である。フーリエ変換の振幅の絶対値の 2 乗を信号 f(x)のパワー スペクトルといい € Φ(ν) =| F(ν) |2 (2.24) で表すことにする。このパワースペクトルは、信号 f(x)にνの周波数成分がパワ ーとしてどれだけの量含まれているかを表している。 図2.3 矩形波(a)のフーリエ変換(b) 以下に、原関数とそのフーリエ変換された関数の間に成り立ついくつかの性 質について示そう。フーリエ変換は数学的な操作であり、これを演算子としてFT で表す。また、フーリエ逆変換を FT-1で表すものとする。また、原関数を英 小文字で、そのフーリエ変換を対応する大文字で表すものとする。フーリエ変 換は、線形操作であるから、次の関係が成り立つ。 € FT[αg(x) +βf (x)] =αG(ν) +βF(ν) (2.25) € FT[g(ax)] = 1 | a |G( ν a) (2.26) α、β、a は、定数である。(2.26)式によると、フーリエ変換においては、元の座 標がa 倍になるとフーリエ変換面における座標は 1/a となることを表している。 次に、座標をa だけシフトした関数のフーリエ変換を考えると € FT[g(x − a)] = exp(−i2πaν)G(ν) (2.27) となる。ここでも、a は定数である。原関数で座標が a だけずれても、そのフー リエ変換の絶対値の形は変わらず、直線的に変化する位相項が付加されるとい うことを表している。このことから、座標移動した関数のフーリエ・パワース ペクトルは、その形は元のままであることを表している。これは、パワースペ クトルの移動不変(shift invariant)と呼ばれている。次に、原関数をフーリエ変換 したものをさらに逆変換すると € FT−1[FT[g(x)]] = g(x) (2.28) となり、元の関数に戻る。ついでながら、フーリエ変換したものを、再度フー リエ変換すると、関数形は元に戻るが、座標が反転する。すなわち € FT[FT[g(x)]] = g(−x) (2.29) となる。2 次元の関数について、関数がそれぞれの変数関数の積で与えられると き、そのフーリエ変換はそれぞれの関数をフーリエ変換したものの積となる。 €
FT[g(x, y)] = FT[gx(x)gy(y)] = FT[gx(x)]FT[gy(y)] (2.30)
この他にも、微分、積分関数のフーリエ変換などに関する有意義な性質など があるが、微分した関数のフーリエ変換について述べておこう。部分積分を使 い、ある関数を微分したものをフーリエ変換すると € df (x) dx −∞ ∞
∫
exp(−i2πνx)dx = i2πνF(ν) (2.31) となる。ここで、f(x)は有限の値の関数であり、さらに € f (±∞) ≈ 0と仮定した。こ の関係を一般化すると € FT[ d n dxn f (x)] = (i2πν) nF(ν) (2.32)が得られる。 2.7 離散フーリエ変換 フーリエ変換は、信号処理、波形処理において非常に応用範囲が広い。しか し、元の変換式に基づいて計算機による数値計算をしようとすると、以下に述 べるように膨大な計算量が必要になる。しかし、幸いなことに、フーリエ変換 を離散的な式として表していくと、同じ計算が随所で出てくる。このため、実 際の計算量を桁違いに節約し、計算の高速化ができる。そのツールとして高速 フーリエ変換(FFT: Fast Fourier Transform)は欠かすことのできないものである。 したがって、本節では、FFT の原理について簡単について述べる。フーリエ変 換自体は(2.20)式で定義されるため、この定義に乗っ取って積分計算を dx の微小 要素について区分数値計算することにより、原関数 f(x)のフーリエ変換である F(ν)の関数形を求めることができる。しかし、積分の指数部分は細かく振動する 成分を含み、さらに無限大積分区間において 2πごとに繰り返す指数部分の値を 毎回計算する必要があり、計算量がきわめて多くなる。特に、デジタル計算で は cos、sin の計算、積の計算には多大な時間を要するため、一般に定義通りに 順番に cos、sin を計算し原関数との積を作り積分計算を行っていたのでは膨大 な計算時間が必要であることがわかる。これらの計算を大幅に軽減するために 考えられた方法が、コンピュータ処理のための高速フーリエ変換である。ここ では、1 次元の関数についての FFT について説明するが、この結果を拡張して 容易に2 次元の場合にも適用できる。 フーリエ変換の式 (2.20)を x、νに対応する微小量Δx、Δνについて離散式で表 すと € F(nΔν) = f (lΔx)exp(−i2πnΔν ⋅ lΔx)Δx l =−∞ ∞
∑
(2.33) のように書き表すことができる。n も整数である。この微小量Δx、Δνについて は任意でよいわけではなく、第 4 章で述べるサンプリングの定理にしたがい、 離散信号が正しく連続な原関数を表現することができるサンプリング間隔にな っている必要がある。ここでは、サンプル間隔はこの定理を満たすものとしよ う。実際には無限和をデータとして作ることはできないため、n、l について同 じN 個の有限個のサンプリング数を考える。Δx=1 とスケーリングし、 (2.33)式 を下記のように表してみる。€ Fn = flWN nl l =0 N −1
∑
(2.34) ただし € Δν = 1 NΔx = 1 N (2.35) € WN nl = exp −i2πnl N ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ = cos 2πnl N ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ − isin 2πnl N ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ (2.36) である。ところで、 (2.34)式において、十分大きな N について € F0 = FN ≈ FN −1と なっている。このことは、fl の関数形にかかわらずフーリエ変換面において (N-1)/2 を境にスペクトルが折り返されていることを表している。すなわち、N/2 以上の成分は負の周波数となっていることを意味している。したがって、実際 に意味のあるのは計算したスペクトルの半分ということになる。 2.8 高速フーリエ変換−FFT 念のため、逆フーリエ変換が離散式としてどのように書けるか導いておこう。 まず、フーリエ逆変換の関係式(2.19)を参考にして € FlWN −nl l =0 N −1∑
がどう計算されるか について考えてみよう。 (2.34)式を代入して € FlWN−nl l =0 N −1∑
= fkWNklWN−nl k =0 N −1∑
l =0 N −1∑
= fk WN(k −n )l l =0 N −1∑
k =0 N −1∑
= fk⋅ Nδkn k =0 N −1∑
= Nfn (2.37) すなわち € fn = 1 N FlWN −nl l =0 N −1∑
(2.38) が得られる。1/N は、フーリエ変換前後でのエネルギー保存に伴う項であるが、 (2.35)式のようなスケーリングをしたことによって出てくる項である。実際の計 算においては、フーリエ変換、逆変換を同じ信号に何回も施す場合にはこの1/N の項を忘れてはならない。なぜならば、この項を忘れると、フーリエ変換の前 後でのエネルギー保存が崩れ、N が大きい場合には丸め誤差などによる計算誤 差が繰り込まれるため正確な計算結果が得られないことがあるからである。 (2.34)式の形は離散フーリエ変換(DFT: Discrete Fourier Transform)と呼ばれる。(2.34)式の計算では、n を 0 から N−1 までとすると、f と W のかけ算を N2回計算 する必要がある。また、(2.36)式において、cos、sin の計算をそれぞれやはり N2 回計算する必要がある。しかし、計算においては同じ値の計算がくり返し行わ れている。たとえば、N=8 を例にとると、F0の計算に € f4W80という項が含まれる が、F2において € f4W88 = f4W8 0、F 4で € f4W816 = f4W8 0、F 6で € f4W824 = f4W8 0のように 同じ項の計算が何回も出てくる。このルールを整理することにより、実際に必 要とされる(2.34)式のかけ算の回数を N2から N 回に減らすことができる。当然 のことながら、N を大きく取ると、計算回数は文字通り桁違いとなり、DFT に 比べ大幅な計算時間の短縮が可能になる。さらに、cos、sin の計算においては、 第1 象限の計算をやり、それらの値に符号をつけることによって他の象限の cos、 sin の値を表現できる。また、cos と sin は独立ではなく、片方を計算すれば他方 は自ずと計算されている。したがって、このcos と sin の計算も実際には 1/8 で すみ、この結果をあらかじめ表(すなわち cos、sin テーブル)にしてコンピュ ータ内に作っておけば、この部分にかかる計算時間を格段に節約できる。 2.9 FFT の実際の計算方法 このような観点に立って、計算回数を節約する実際の計算方法である FFT と して、いくつかの方法が提案されている。いずれの方法においても、基本とな るのは図 2.4(a)に示すバタフライ演算である。ここでは、クーリー・ターキー (Cooley-Tukey)アルゴリズムにそった FFT について説明してみよう。図 2.4(b)は N=8 とした例である。2 入力 2 出力のバタフライ演算を基本とすると、N=23=8 の場合、3 段階のバタフライが必要となる。2 つの入出力のみについて考えると、 図2.4(a)に示すように、A からそれぞれへの出力は重みのない加算として取り扱 われる。しかし、B からのそれぞれへの出力においては、矢印の下に書かれた フーリエ変換の回転子Ws、Wtがそれぞれ重みとしてかけられ、その出力への加 算となる。これを基本として、N=8 に拡張したのが図 2.4(b)である。たとえば、 Wmはm=8 ごとにくり返す値を持ち、Wmでm の次数についての等価な W の値を 考慮すると、この図のアルゴリズムからF3の計算は € F3 = f0+ f1W 3 + f2W 6 + f3W 6+3 + f4W 4 + f5W 4 +3 + f6W 4 +6 + f7W 4 +6+3 = f0W 0 + f1W 3 + f2W 6 + f3W 9 + f4W 12 + f5W 15 + f6W 18 + f7W 21 (2.39)
となり、F3について正しく(2.34)式が計算できていることがわかる。バタフライ 演算からわかるように、FFT の基数は 2 であり、したがって、計算すべきデー タの数N としては通常 2n個(n 整数)とする。図 2.4 (b)に示すように、データ 処理の矢印のフローとしては上下中央部を境にして対称な形であるが、入力デ ータf は上から順番に整列しているのに対し、出力データ F をみると偶数と奇数 図 2.4 FFT におけるバタフライ演算(a)と計算のフロー(b)
番目のデータが上下にわかれており、さらにそれぞれで順番が入れ替わってい る。このため、実際の FFT プログラムにおいては、正しい出力波形を順番に表 示するために並べ替え(これをスワップという)のプロセスが必要になる。2 次 元画像などに対する FFT でも、ここでの議論を拡張すればよいが、どの計算方 式をとっても、整列した入力データについて計算された変換データの並び順は あるルールで入れ替わっており、出力データのスワップも重要なソフトウエア の部分を占めることに注意しておこう。FFT は、コンピュータを用いた周波数 解析においてなくてはならない道具である。 課題 2 フーリエ変換の関係式について、(2.26)、(2.27)、(2.28)、(2.29)、(2.31) 式が成り立つことを確かめよ。
3.線形システムと伝達関数 3.1 線形システム 時間変化する信号は、線形システムを伝搬して伝わることが多い。ここでは、 そのような信号の伝搬について入力と出力の関係について調べてみよう。多く の物理系、工学系では、入力信号、信号伝搬システム、出力信号の三つが重要 であり、このうち二つの性質がわかかっているものとして、残りの一つを求め るというのが通常の問題設定になっている場合が多い。たとえば、信号が伝搬 するシステムの性質がわかかっており、出力信号を得たとき入力信号がどのよ うな信号であるかを求めるいわゆるセンシングシステムの場合などがそれであ る。また、入力が分かっていてその出力を得たときには、システム関数がどの ようであったかを知ることができる。一方、入力も分かっていて、システム関 数が分かっている場合に、出力がどうであるかを知ることも重要である。この 場合は入力、システムが分かっているわけであるから、素直にそのまま出力を 得ればよいことになるが、大規模システムなどの場合には、システムを設計し たときに出力どのようになるかの事前シミュレーションが必須になる。これが この場合に当たる。 図3.1 のような線形システムで、入力信号を f(x)、出力信号を g(x)として、入 出力の間の関係を求めてみよう。いま、入力信号をデルタ関数を使って € f (z) = f (ξ)δ(z −ξ)dξ −∞ ∞
∫
(3.1) と置き換える。(3.1)式のように、ある関数(ここではデルタ関数)の座標をシフト したものと他の異なる関数との積を積分した形は、畳み込み積分(convolution)と 呼ばれる。積分変数ξについて、相関のところで述べる(4.1)式の積分計算におけ る関数f の変数はξ−x で表されているのに対し、(3.1)式ではδ関数の変数は z−ξと なっており、変数の符号が反転していることに注意しよう。 3.2 線形システムにおけるたたみ込み積分 システム応答を数学的な演算子としてS で表してみよう。 (3.1)式に対して S の演算をほどこすと、出力関数g(x)は € g(x) = S[ f (z)] = f (ξ)S[δ(z −ξ)]dξ -∞ ∞∫
(3.2)となる。ここで、出力の座標は元のz と同じである保証はなく、時間遅れなどを 含む可能性があり、一般的にz とは異なる座標 x とした。また、S は座標 z に対 する演算子であり、積分の順序の入れ替えをしている。(3.2)式のデルタ関数に 対するS の演算を € S[δ(z − ξ)] = h(x,ξ) (3.3) と置いてみよう。関数h(x)は、デルタ関数に対する線形システムの応答関数であ り、インパルス応答と呼ばれる。すなわち、関数h(x)は、時間変化する現象の例 では短時間のパルスを入力としたとき、システムを伝搬した後の信号がどのよ 図3.1 入出力線形システム うに広がるかを表している。理想的なシステムでは、インパルスの入力に対し て出力として入力と同じ波形のインパルスが得られる。このとき、システムの 応答関数はデルタ関数で表されることになる。しかし、一般的な入出力におい ては、パルス入力は出力側でより広がった幅をもつことになる。システムのイ ンパルス応答がわかかったとして、(3.3)式を使うと、一般的な入力関数 f(ξ)に対 して € g(x) = f (ξ)h(x,ξ)dξ −∞ ∞
∫
(3.4) として出力関数が得られる。 線形システムの多くは、時間あるいは空間に対して不変な系となっており、 このような系ではシステム関数は € h(x,ξ) = h(x −ξ) (3.5) のように座標の差だけの関数として与えられる。このとき、信号出力は€ g(x) = f (ξ)h(x −ξ)dξ = f (x) ∗ h(x) −∞ ∞
∫
(3.6) のように、入力信号とシステム応答の畳み込み積分で与えられることになる。 関数f(x)と h(x)の畳み込み積分は、(3.6)式の最後の式の記号*を使ってその演算 を表すことがある。ここでも、畳み込み積分を同様の記号で表すことにする。 すでに述べたように、(3.6)式は(4.1)式の相互相関の式と似た形になっているが、 積分の中の二番目の関数中で、座標の差のとり方が異なっていることに注意し たい。相関のときには、関数 f(ξ)とξについて座標を x だけずらした関数 g(ξ-x) との積を取り、これを積分したものである。これに対し畳み込み積分では、関 数f(ξ)とξについて座標を x だけずらし裏返しにした関数 h(x-ξ)との積を取り、こ れを積分したものである。 3.3 たたみ込み積分の意味 ここで、 (3.6)式の意味を考えてみたい。この式では、積分の中の関数 h は h(ξ-x) ではなくh(x-ξ)となっている。関数形が h(ξ-x)の場合には、たたみ込み積分では なく、後で述べる相互相関関数になる。どちらでもよさそうに思えるが、物理 的な意味からいうと、後者でなければならない。図3.2 を見ていただきたい。わ かりやすくするために、最初に離散的な関数で説明しよう。n=0 以前には信号は ないものとする。n=2 の時刻における(3.6)式の出力に相当する g がどうなるかを 図 3.2 離散関数 f と h のたたみ込み 考える。まず、f(0)の信号は n=2 において図3.3 たたみ込みの概念 (a)入力関数、(b)システム関数、(c)反転システム関数 € g0(2) = h(m) f (0) m =0 2
∑
(3.7) である。同様に、f(1)、f(2)の信号について € g1(2) = h(m −1) f (1) m =0 2∑
(3.8a) € g2(2) = h(m − 2) f (2) m =0 2∑
(3.8b) となる。h(m)で m が負となるときには h(m)=0 であることに注意しよう。したが って、n=2 において、全出力 g(2)は € g(2) = h(n − m) f (m) m =0 2∑
(3.9) と与えられる。この関係を連続関数として表すと、図3.3 に示すように、図 3.3(b)のレスポンス関数を垂直軸に対して図 3.3(c)のように反転させ、図 3.3(a)の入力 関数に座標を移動させながら積を取り、重なり合う領域の積分計算を行うこと になる。すなわち、(3.6)式になる。これがたたみ込み積分の中身である。相関 関数とは類似の演算ではあるが、物理的意味合いは異なる。 3.4 伝達関数 入力信号f(x)を電圧などの時間変化する関数にたとえれば、それに付随する信 号の周波数成分を考えることが出来る。この入力信号のフーリエ変換を F(ν)で 表すことにしよう。そうすると、図3.1 において、システム関数 h(x)、出力関数 g(x)についても、それぞれのフーリエ変換 H(ν)、G(ν)を定義することが出来る。 (3.6)式に示したように、入力、システム、出力の関数はたたみ込み積分によっ て関係づけられることがわかった。それでは、フーリエ変換された周波数面で 定義されるそれぞれの関数同士の関係はどのようになるであろうか。これらを 調べることにより、単に実関数面と周波数面での対応関係を明らかにするにと どまらず、場合によっては周波数面で信号解析を行うことにより、実関数面で 処理を行うよりも信号処理、解析を容易にできるなどの利点が明らかにされる。 フーリエ変換は線形操作であるから、実関数面、フーリエ面、どちらの面での 信号解析を行うことができる。それぞれの結果をフーリエ変換あるいは逆変換 することによって、必要とされる時間情報、周波数情報を最終的に得ることが 可能になる。 (3.6)式でシステム関数 h(x)がわかかっているとき、出力 g(x)から入力関数 f(x) を求めるためには、この積分方程式を解く必要がある。これは、一般的に手間 のかかる方法であるが、それぞれの関数のフーリエスペクトルを使うと、積分 方程式を解かずに入力関数f(x)を求めることができる。出力関数 g(x)のフーリエ 変換は € G(ν) = { f (ξ)h(x −ξ)dξ}exp(−i2πνx)dx −∞ ∞
∫
−∞ ∞∫
= f (ξ)h(x −ξ)exp{−i2πν(x −ξ) − i2πνξ}dξdx −∞ ∞∫
−∞ ∞∫
= [h(x −ξ)exp{−i2πν(x −ξ)}dx] f (ξ)exp(−i2πνξ)dξ −∞ ∞∫
−∞ ∞∫
= H(ν)F(ν) (3.10) となる。関数の畳み込み積分は、対応するフーリエ面においてはそれぞれの関 数をフーリエ変換したものの積で表されることがわかかる。このことから、システム応答のフーリエ変換H(ν)がわかっており、出力のフーリエ変換 G(ν)が得 られたとき、入力関数のフーリエ変換はF(ν)=G(ν)/H(ν)によって計算される。求 められた F(ν)を逆フーリエ変換すると、元の入力関数 f(x)が得られる。ただし、 H(ν)にはかならず零点が含まれ、零点となるνの値においては F(ν)を正確に求め ることができない。このとき、信号回復として用いる € T(ν) = 1/ H (ν)は逆フィル タと呼ばれる。元信号を周波数面において正しく再生する問題については、5.2 節において詳しく学ぶ。H(ν)をシステムの伝達関数(transfer function)という。ま た 、 入 力 信 号 の 周 波 数 面 に お け る パ ワ ー は 、 入 力 パ ワ ー ス ペ ク ト ル € ΦF(ν) =| F(ν) |2として定義される。出力信号のパワースペクトルを € ΦG(ν) =| G(ν) |2 とすると、これらのパワースペクトル間について、(3.10)式より入出力のパワー スペクトルの関係は € ΦG(ν) =| H(ν) |2ΦF(ν) (3.11) となる。 課題3 € H(ν)F(ν)を逆フーリエ変換して、結果が € h(x) ∗ f (x)となることを証明せ よ。
4.信号の相関、パワースペクトル、サンプリングの定理 4.1 相関関数 関数 f(x)の自己相関関数は、ある x での値 f(x)が異なる x'における値 f(x')とど の程度まで関係があるか、あるいはどの程度似ているかを表す関数である。相 関関数は、関数のフーリエ変換とも深い関係があり、信号の性質を評価すると きにしばしば用いられる。また、関数f(x)とこれと異なる関数 g(x)との間の関係 を表すものとして、相互相関関数があり、これはお互いの関数がどの程度似通 っているかを評価するために用いられる。これについても後でのべる。ここで は、先に自己相関関数の定義、性質について調べる。相関関数は、前節で述べ たシステム応答を表す畳み込み積分とも深い関係がある。 自己相関関数は € R(x) = lim T →∞ 1 T f (ξ) f *(ξ − x)dξ −T / 2 T / 2
∫
(4.1) で定義される。T は適当にとられた座標 ξ の積分間隔であるが、以下では相関積 分は有限値に収束するものとし、(4.1)式の極限操作はいちいち書かずに積分は 無限区間とし、簡略化して表すこととする。また、相関関数は、積分の代わり に6 節で述べる集合平均の記号<・>を使い € R(x) =< f (ξ) f*(ξ − x) > (4.2) と表されることもある。集合平均と(4.1)式の平均(x を時間に見立てれば時間平 均)とは意味が異なるが、たいていの物理的信号の場合には統計関数に関して エルゴード性が成り立ち、6 節でも述べるように集合平均と積分平均は同じ値を とることが示される。関数f(x)は、一般的に複素関数として定義しているが、f(x) が実関数のときには € R(x) = f (ξ) f*(ξ − x)dξ = f (ξ + x) f*(ξ)dξ −∞ ∞∫
−∞ ∞∫
= R(−x) (4.3) となり、偶関数となる。相関関数として(4.2)式のように書いたときには、(4.1) 式の意味であるとする。また、x=0 における相関関数の値は € R(0) = | f (ξ) |2 dξ −∞ ∞∫
(4.3) であり、R(0)は、信号の全エネルギーを表している。4.2 相関関数の例 以下に、いくつかの相関関数の例をあげる。図4.1(a)に示す矩形パルス波の相 関関数は € R(x) = rect(ξ)rect(ξ − x)dξ = −∞ ∞
∫
−1/ 2−xdξ =1− x x ≥ 0 1/ 2∫
dξ =1+ x x < 0 −1/ 2 1/ 2+x∫
⎧ ⎨ ⎪ ⎩ ⎪ (4.4) と計算され、図4.1(b)のように表される。もう一つの例として、図 4.2(a)に示し たようなランダム関数の相関に関して述べておこう。ランダム関数は、解析的 な形で関数を表すことはできないが、(4.1)式の相関関数の定義式に従って、相 関関数を数値的に計算することができる。ランダム雑音の場合、図4.2(b)の例よ うに、一般にガウス型に近い分布の相関関数が得られる。通常の信号処理で現 れるランダム関数といえども、ランダムと見なせる時間について有限な時間幅 を持っており、時間が経つと十分ランダムではあるが、短い時間では相関を持 つ信号であることが多い。この関数の幅は、ランダム現象の平均的な帯域を表 しており、たとえばこれからガウスゆらぎをするブラウン運動粒子の固有時間 変化などを調べることができる。 図4.1 矩形波の相関 (a)矩形波 (b)相関関数 図4.2 ランダム関数の相関 (a)ランダム関数 (b)相関関数 相互相関関数も、自己相関関数と同様に定義でき、異なる二つの関数f(x)と g(x) の相互相関関数は € R(x) = f (ξ)g*(ξ − x)dξ −∞ ∞∫
(4.5)と表される。相互相関関数は、お互いの関数が座標 x に関してどの程度似よっ ているかを表している。 4.3 相関関数とパワースペクトル フーリエ変換と相関関数との間の関係について述べておこう。(2.24)式のパワ ースペクトルの定義を使い、パワースペクトルの逆フーリエ変換を計算すると € Φ(ν)exp(i2πνx)dν = f (x)exp(−i2πνx)dx −∞ ∞
∫
−∞ ∞∫
−∞ ∞∫
2exp(i2πνx)dν = { f (x1)exp(−i2πνx1) f * (x2)exp(i2πνx2) −∞ ∞∫
−∞ ∞∫
−∞ ∞∫
dx1dx2}exp(i2πνx)dν = f (x1) f * (x2)exp{i2πν(x2− x1+ x) −∞ ∞∫
−∞ ∞∫
−∞ ∞∫
}dνdx1dx2 = f (x1) f *(x 2)δ{x2− (x1− x)}dx1dx2 −∞ ∞∫
−∞ ∞∫
= f (x1) f *(x 1− x)dx1 −∞ ∞∫
= R(x) (4.6) となる。すなわち、パワースペクトルのフーリエ変換は、相関関数となること がわかる。したがって、パワースペクトルと相関関数は、同じ情報を含むこと になる。これを、ウイーナー・キンチンの定理(Wiener-Khintchine theorem)とい う。また、この式でx=0 と置くと € Φ(ν)dν = | F(ν) |2 dν = | f (x) |2 dx −∞ ∞∫
−∞ ∞∫
−∞ ∞∫
(4.7) が得られ、フーリエ変換面におけるエネルギーと実面におけるエネルギーは同 じ、すなわちどのような変換面においても、エネルギーの保存が成り立つこと を示している。これを、パーサバルの定理(Parseval theorem)という。ここでは証 明はしないが、一般化されたパーサバルの定理について述べておく。関数 f(x) とg(x)について、一般化パーサバルの定理は € f (t)g*(t)dt −∞ ∞∫
= F(t)G*(t)dt −∞ ∞∫
(4.8) のように表される。関数 f(x)、g(x)は複素関数に拡張されており、大文字はそれ ぞれのフーリエ変換を表す。この式において、f(x)=g(x)とすると(4.7)式に等しく なることは言うまでもない。 4.4 信号のサンプリング たとえば、時間的に変化する物理信号は、時間に対して切れ目なく連続した信号となっている。このような信号を、計算機などを用いて解析しようとする と、時間について離散的な間隔でこの信号を読み取らなければならない。では、 いったいどのくらいの時間間隔で信号をサンプルすればよいのだろうか。時間 に対して細かく信号をサンプルすれば、より元の信号を忠実に再生できるであ ろうことは間違いない。しかし、データ量は膨大になる。これに対してあまり サンプル間隔が開き過ぎると、今度は元の情報が失われてしまう。より、少な いデータで元の信号についての最大(数学的には完全な)の情報を与えようと するのが、ここで述べるサンプリングの定理である。このことは、空間的な画 像についても同様な議論として成り立つ。すなわち、連続する空間画像につい てもやはり最適なサンプリング画素の単位が存在する。 図4.3 関数のサンプリング サンプリングの定理は、信号の最大帯域幅がわかっているとき、この帯域幅 に関係した一定の時間間隔で信号をサンプリングすれば、元の信号を正しく再 生することができるというものである。サンプリングの定理は、フーリエスペ クトルと関係が深く、フーリエ変換の性質を使うことにより説明ができる。以 下では、この定理を証明する前に関数のサンプリングについて述べよう。信号 を関数f(x)で表し、図 4.3 に示すように、この信号を x について等間隔 X でサン プリングする場合を考えてみよう。サンプリングされた関数を fs(x)とすると、 fs(x)は € fs(x) = comb( x X) f (x) (4.9) と書ける。ここで、comb(x)はコム関数(櫛の歯関数)と呼ばれ € comb(x) = δ(x − n) n=−∞ ∞
∑
(4.10) のように定義される。実際の信号のサンプリングとしては、時間軸だけでなく、 たとえば電気回路における電圧信号であれば、電圧値についても AD 変換によって離散的な値としてサンプルされる。しかし、ここでは電圧値は忠実なアナ ログ値として取得されるものとし、時間軸だけについて信号がサンプルされる 場合を考える。 4.5 サンプル関数のフーリエ変換 サンプルされた信号は、そのフーリエ変換面でどのように表されるかについ て考えてみよう。(4.9)式から、サンプリングされた関数 fs(x)のフーリエスペクト ルは € Fs(ν) = FT[comb( x X)]∗ F(ν) (4.11) と書ける。comb 関数のフーリエ変換は、やはり comb 関数になることが、次の ようにして容易に確かめられる。(4.10)式のフーリエ変換は € comb(x)exp(−i2πνx)dx = −∞ ∞
∫
exp(−i2πnν) n=−∞ ∞∑
(4.12) となる。ここで、くり返しの関数についてフーリエ級数展開のところで行った ように、comb 関数を級数展開して € comb(x) = anexp(−i2πnx) n=−∞ ∞∑
(4.13) のように書く。ここで、係数anは € an= comb(x)exp(−i2πnx)dx = 1 −1/ 2 1/ 2∫
(4.14) となり、(4.10)の右辺の関数形は(4.13)式と等しくなる。このことを使うと、(4.11) 式は € Fs(ν) = δ(ν − n X) ∗ F(ν) = F(ν − n X n=−∞ ∞∑
n=−∞ ∞∑
) (4.15) と計算される。ここで、(4.11)式の comb 関数のフーリエ変換は、フーリエ変換 の性質を使うと € FT[comb(x X)] = Xcomb(Xν) = Xδ(Xν − n) (4.16) となるが、(4.15)式を導くために € δ(Xν − n)g(ν )dν = 1 X δ(ν − n X −∞ ∞∫
−∞ ∞∫
)g(ν)dν (4.17)となる関係を用いた。すなわち、先にも述べたように、δ 関数は積分の中におい て意味を持つ関数であるということを忘れてはならない。 図4.4 サンプリング関数のスペクトル 4.6 サンプリングの定理 さて、サンプルされた関数のフーリエ変換はわかったので、この関数の中に どのように元関数の情報が含まれているかについて考えよう。(4.15)式を見ると、 サンプリングされたスペクトルは、図4.4 に示すように、元の信号のスペクトル F(ν)をν=1/X の間隔で並べたものであることがわかる。はじめにも述べたように、 信号はある最大の帯域を持っているものとしているので、この帯域B が 1/2X の 幅よりも小さいときには、サンプリングされたスペクトル面で、互いに隣のス ペクトルと重なることはなく、ν=0 の周りにおいて元の信号のスペクトルのみを 分離することが可能になる。したがって、分離できる上限のサンプリング間隔 の条件は € X ≤ 1 2B (4.18) である。すなわち、元信号のサンプリング間隔をX=1/2B とするとき、最もよい 効率でデータのサンプリングができることがわかる。 さて、図4.4 に示したようなサンプリングされた信号のスペクトルから、元の 信号のスペクトルのみを取り出すには、次のようなフーリエ面でのフィルタを 通すとよい。 € H(ν) = rect( ν 2B) (4.19) X=1/2B のとき、(4.19)式で表されるフィルタ関数をサンプルされた関数のフーリ エ変換式(4.15)にかけると € Fs(ν)H (ν) = F(ν) (4.20)
となり、正確に元の信号のフーリエ変換が得られる。 (4.20)式をフーリエ逆変換すると、右辺は元の信号 f(x)となるが、左辺はどのよ うになるであろうか。(4.20)式右辺を逆フーリエ変換すると € f (x) = [comb( x X) f (x)]∗ h(x) = { δ( x X − n) f (x)} ∗{2Bsinc(2Bx)} n=−∞ ∞
∑
= f ( n 2B)sinc{2B(x − n 2B)} n=−∞ ∞∑
(4.21) となる。(2.64)式によると、図 4.5 に示すように、X=1/2B の間隔でサンプリング した信号に、サンプリングの間隔に等しい広がりを持つ sinc 関数を掛けて足し 合わせることにより、元の信号が正しく再生されるということがわかる。これ を、シャノン(Shannon)のサンプリングの定理と呼ぶ。 図4.5 サンプリング関数からの元信号の復元 4.7 サンプリングの定理の意味 実際のシステムでは、伝送できる周波数帯域、計測などにおいては検出可能 な周波数帯域が限定されているため、サンプリングの定理の適用は非常に有効 である。図4.5 からわかるように、信号帯域 B に対し € X << 1/2Bとなるサンプリ ングをすると、(4.19)式のフィルタ関数の幅に対しサンプルされた関数幅がかな り小さくなるため、無駄な周波数情報を使って原信号を再生することになる。 すなわち、必要とされる情報以上に多くのデータを取得することになる。これ をオーバーサンプリングというが、効率的な信号サンプリングとしては望まし くない。一方、信号帯域B に対し € X >>1/2Bとなるサンプリングをすると、サン プルされた関数のフーリエ変換の幅が式(2.62)のフィルタ関数の幅を超え隣り 合うスペクトルと重なるため、フィルタ関数をかけて信号再生しても、原関数 を正しく再生することはできない。この場合は、信号のサンプルが粗すぎて、 元の情報の欠落が起こる場合である。これをアンダーサンプリングというが、これも信号サンプリングとしては望ましくない。したがって、先に示したよう に € X =1/2Bとなるように信号サンプリングをするのが最適なサンプリング条件 であると言える。しかし、現実の信号には一般に雑音が含まれるため、実際の 信号サンプリングにおいて € X =1/2Bでのサンプリングを行っても正しく信号再 生をすることは難しい。雑音がある場合にどのような信号サンプリングが望ま しいかは、雑音のレベルによって条件は変わるため、一般的な議論は難しい。 雑音がある場合には、理想的なサンプリング条件から求められる € X =1/2Bより も多少サンプリング間隔を小さくして、多めの信号サンプリングを行う必要が ある。 課題4 一般化パーサバルの定理(4.8)式を証明せよ。