3. フーリエ変換 3. 周期的な複雑な波形 (t) si(ωt), (t) si(ωt), (t) si(3ωt) のグラフを図 3 に示す 単純にこれらの波形を重ね合わ せると (t) si(ωt) + si(ωt) + si(3ωt) は右図のように複雑な波形となる この合成波の時間方向の移 動は見られない ( 時間方向を波の位相と呼ぶ ) しかし 振幅の変調が見られる 3 3Hz (t) Hz (t) - Hz -..5..5. ime (sec) -3..5..5. ime (sec) 図 3 周期関数 一方 (t) cos(ωt) + si(ωt) は 図 3 のようになる 振幅の変調はないが 位相のズレが見られる 例えば ( t) A cos( ωt) + B si( ωt) A B A + B cos( ωt) + si( ωt) A + B A + B (t) - A + B A + B si t { si( φ) cos( ωt) + cos( φ) si( ωt) } ( ω +φ) -..5..5. ime (sec) (3 ) 図 3 si(ωt) + cos(ωt) となるので 位相 φだけずれる つまり si と cos の適当な組み合わせで任意の周期関数を表すことができる この考えに基づいて式に表すと (t) a + a cos(ωt) + a cos(ωt) + a3 cos(3ωt) + + b si(ωt) + b si(ωt) + b3 si(3ωt) + (3 ) となる ( フーリエ級数 )
3. 波形合成の意味 周期的な波形は振幅 ( a, b ) と位相 ( ω π ) で特徴付けられることがわかったが 実際にはどういう意味なのか? 例えば 色は 3 原色 ( 赤 (R) 緑(G) 青(B)) の混ぜ方ですべての色を表すことができる 色は周波数 ( 波長 ) で決まり 色の明るさ ( 強度 ) は振幅で決まる つまり a, b, ω の組み合わせであらゆる色が作り出される つまり 料理では砂糖 塩 醤油を何グラムずつ配合するかによっていろいろな味ができるのと同じである 音は 基本振動数 倍振動数 3 倍振動数 の割合で波形が決まり いろいろな音の 音色 が生み出される 3.3 波の構成要素 (3 ) 式のように 簡単な波の重ね合わせで複雑な波形を作れることがわかったが ここでは逆に複雑 な波形から 基本的な波の構成要素 (ω, ω, に対応する a, a, b, a, b, ) を求める あ る音の複雑な繰り返しパターンの基本振動数を ( ω π ) とする 図 3 の場合は 周期 (sec) なので 周波数 (Hz) となる まず si, cos の積分を計算してみる ( ω t) si cos( ω t) dt ( ω π ) ω si ( ω t) cos dt ω ( ω t) となるので cos( ω t) dt, ( 3ω t) dt,, ( ω t) dt, ( 3ω t) dt cos が求まる 結局 (3 ) 式の積分は ( t) dt a si になり si, a ( t) dt (3 3) ここで 複雑な波形の 周期分を積分すれば a が求まることがわかった 次に a, b, a, b, を 計算する その前に { a cos( ω t) a cos( ω t) } dt の計算をしてみる cos どうしの掛け算の積分 ( ) ( A+ B) + cos( A B) cos cosa cosb の公式を用いて aa { a cos( ω t) a cos( ω t) } dt cos( 3ω t) dt+ cos( ω t) dt
cos どうしの掛け算の積分 ( ) ( A) cos + cos A から a となる a dt + ( ω t) cos( ω t) dt cos( ω t) dt cos a 3 si どうしの掛け算の積分 ( ) cosa+ B 同様に sia sib b si ( ) cos( A B) bb を利用して ( ω t) b si( ω t) dt cos( 3ω t) dt cos( ω t) dt 4 si どうしの掛け算の積分 ( ) ( A) cos si A から となる b b si( ωt) si( ωt) dt dt cos( ωt) dt b 5 si cos の掛け算の積分 ( ) si sia cosb ( A+ B) + si( A B) ba b si ( ω t) a cos( ω t) dt si( 3ω t) dt+ si( ω t) dt 6 si cos の掛け算の積分 ( ) si A cosa si( A) ba b si ( ω t) a cos( ω t) dt si( ω t) dt まとめると cos どうしの掛け算の積分 ( ) と si どうしの掛け算の積分 ( ) の場合だけ 積分値 が にならない b siωt b siωt a cosωt a cosωt siωt b / cosωt a /
これらのことを考慮して (3 ) 式の両辺に cosωt をかけて積分すると ( t) cos( ω t) dt a cos( ω t) dt+ a cos( ω t) cos( ω t) dt+ a cos( ω t) cos( ω t) L dt+ + ( ω t) cos( ω t) dt+ b si( ω t) cos( ω t) L b si dt+ a a ( t) cos( ωt) dt 一般に 次のような関係が求められる a ( t) cos( ωt) dt (3 4) b ( t) si( ωt) dt (3 5) 結局 (3 3) 式 (3 4) 式 (3 5) 式によって フーリエ級数の係数の求め方がわかった 3.4 数値積分 測定データに対して 積分をしないとフーリエ係数が求まらないので まずは 数値積分のプログラムを作る 実際の測定データは離散データになるので 積分から和になる ( x) dx ( ) x x (3 6) si θ..5 Δx が小さくなれば連続関数の積分値に近づいていく 実際に siθ ( θ π / ) の離散データ 数値積分を行ってみる 図 3 3 のように 通りの近似 の仕方がある さらに近似度を上げるには 棒状近似 から台形近似 ( シンプソンの台形近似 ) にする ( x ) + ( x ) ( x) x x (3 7) 3 つの場合の計算プログラムは次のようになる 変数 Sum, Sum, Sum に計算結果が代入される シンプ ソンの台形公式 (Sum) が真の値に近いことがわかる また dx を小さくしていくと Sum, Sum, Sum の差が 小さくなっていく 離散データの 3 点を使って 次関数近似 計算もあるが 少し複雑になるので ここでは述べない si θ. 3 4 5 6 7 8 9..5 θ (deg.). 3 4 5 6 7 8 9 θ (deg.) 図 3 3 siθ の数値積分 3
Dim pi As Double Sub iteg() pi 4# * At(#) ' π 円周率 X_max 9 ' 積分範囲 ( ) dx ' X #: X dx I : Sum #: Sum #: Sum # Do I I + Y Si(X * pi / 8#) Y Si(X * pi / 8#) Sum Sum + dx * Y Sum Sum + dx * Y Sum Sum + dx * (Y + Y) / # Cells(I, ) X Cells(I, ) Y Cells(I, 3) X Cells(I, 4) Y X X + dx: X X + dx Loop While X < X_max Cells(I +, ) X_max / dx Cells(I +, ) Sum * pi / 8# Cells(I +, 3) Sum * pi / 8# Cells(I +, 4) Sum * pi / 8# Ed Sub 3.5 離散フーリエ変換 (DF: Discrete Fourier rasorm) (3 ) 式を書き直すと ( t) a + { a cos( ωt) + b si( ω )} t (3 8) i θ e θ iθ + e e オイラーの公式より e i cosθ+ isiθ cosθ siθ ( t) a + a iωt iωt iωt iωt ( e + e ) b ( e e ) + i e i i θ iθ を (3 8) 式に代入して ( t) a + iωt ( a ib ) e ( a ib ) + + e iωt c a ib c a + ib とすると iωt ce iωt iωt iωt ( t) c + { ce + c e } c + ce + ( ) t c e iωt 4 (3 9)
また c a ib (3 4) 式と (3 5) 式より c ( t) { cos( ωt) i si( ωt) } ( t) e iωt dt dt (3 ) プログラムで Hz の複雑な波形を作り フーリエ変換してみる (3 8) 式の が4までの級数でそれぞれの係数の値を以下のようにすると 図 3.4 のようになる a.5 a.5 b.8 a. b -.4 a3 -.7 b3. a4 -. b4.3 (t) 4 3 - - -3..5. t (sec.) (t).5 +.5 cos(ωt) +. cos(ωt) -.7 cos(3ωt) -. cos(4ωt) +.8 si(ωt) -.4 si(ωt) +. si(3ωt) +.3 si(4ωt) 図 3 4 複雑な波形 (Hz) 計算結果は次のようになる 初期値 N a b c.5.5.5.8.47699. -.4.367 3 -.7..353553 4 -..3.68466 計算値 N a b c.486.486.36 -.4.46443.86..776 3 -.364 -.5.36748 4 -.64 -.5.6357 5
Dim pi As Double, D As Iteger, DF As Iteger,,, X(5) As Double, Y(5) As Double Dim a(, ) As Sigle Sub iit_wave() pi 4# * At(#) ' π 円周率 ' Frequecy (Hz) # / ' 周期 (sec) dt / 5 ' t t_max ' 積分範囲 ( 周期 ) ~t~ DF 4 ' a.5 ' 初期値 a(, ).5: a(, ).: a(, 3) -.7: a(, 4) -. a(, ).8: a(, ) -.4: a(, 3).: a(, 4).3 '------------------------------------------------------------- t #: I Do I I + X(I) t Y(I) a For K o DF Y(I) Y(I) + a(, K) * Cos(K * # * pi * * t) + a(, K) * Si(K * # * pi * * t) Next K Cells(I, ) X(I) Cells(I, ) Y(I) t t + dt Loop While t < t_max + dt D I Cells(, 4) Cells(, 7) a For I o DF Cells(I +, 4) I Cells(I +, 5) a(, I) Cells(I +, 6) a(, I) Cells(I +, 7) Sqr(a(, I) ^ + a(, I) ^ ) / # Next I '**************** DF '**************** Ed Sub Sub DF() dt X() - X() For N o DF c #: s # For I o D c c + Y(I) * dt * Cos(-N * # * pi * * X(I)) s s + Y(I) * dt * Si(-N * # * pi * * X(I)) Next I c c / s s / Cells(N + 8, 4) N Cells(N + 8, 5) c Cells(N + 8, 6) s Cells(N + 8, 7) Sqr(c ^ + s ^ ) Next N Ed Sub 6
3.6 周期のわからないフーリエ変換 これまでは周期のある波形のフーリエ変換であったが 実際の波形は周期がはっきりしていないものが多い ( 図 3 5) それならば周期を無限大と考えれば 周期があいまいな複雑な波形の周期を決める必要がなくなる (3 ) 式は 周期分の積分だから次式のように書き換えられる Curret (µa) 5 ( t) / iωt c e dt (3 ) / この式を (3 9) 式に代入すると 6 8 4 3 36 time (sec.) ω ω ( ) / i t i t t ( t) e dt e 図 3 5 実際の波形 / とすると / Δ となる の連続関数 (ω ω) になるので和 ( Δ) は積分 (d) に なる ( t) lim / ( t) iωt iωt πit { e dt} e { ( t) e dt} πit e d / ( ) ( t) G πit e dt とすると (3 ) ( t) { G ( ) } πit e d (3 3) (3 ) 式をフーリエ変換 (3 3) 式を逆フーリエ変換 と呼ぶ 実際に簡単な関数 ( 図 3 6) のフーリエ変換を計算し てみる (t) (t) ( -Δt/ t Δt/ ) (t) ( t < -Δt/, Δt/ < t ) この (t) を (3 ) 式に代入すると ( t < -Δt/, Δt/ < t ) の範囲で になるので ( -Δt/ t Δt/ ) の範囲だけを図 3 6 考えればよい ( ) t / πit πit π e dt e { e i t π e i t } G t/ πi t / t / πi -Δt/ Δt/ t 7
e siθ G e i i θ iθ siπ t π なので 結局 siπ t π t ( ) t siθ となる また lim θ θ の関係を使うと のとき G( ) t t となる 実は Δt に したとき G() となるデルタ関数である つまり 一定関数のフーリエ変換はデルタ関数となる G() プログラムを作って確かめてみる 5 4 3 図 3 7 - -5-4 -3 - - 3 4 5 (Hz) t t t4 3.7 周期がわからない離散フーリエ変換 (DF: Discrete Fourier rasorm) Dim pi As Double Sub iit_wave() pi 4# * At(#) ' π 円周率 '-------------------------------- dt ' t (sec) '-------------------------------- div dt / 5 ' _max 5 ' (Hz) d _max / 5 -_max I Do I I + t -dt / c : s Do c c + div * Cos(-# * pi * * t) s s + div * Si(-# * pi * * t) t t + div Loop While t < dt / + div Cells(I, ) Cells(I, ) c + d Loop While < _max + d Ed Sub (3 ) 式の積分でフーリエ変換の計算が行われるが プログラムを使っての数値計算では図 3 8 のような短冊の和で計算する 短冊の幅を τとすると t kτとなる (kは整数) G πit ( ) ( t) e dt ( kτ) N k e πikτ 短冊が合計 N 本あるとすると Nτとなる 基本周波数は /なので 倍周波数は τ (t) 5 4 3..5. t (sec.) 8
N τ N ( ) π N G τ ( τ) ik/ G k e (3 4) Nτ k 実際に図 3 4 の基本周波数 (Hz) の 4 の計算をプログラムで計算すると.3 図 3 9 のようになる G() は,, 3, 4 で値を持ち (t) の c, c, c3, c4 の値に比例している G().5..5 3 4..5. 3 4 5 6 (Hz) 図 3 9 9
Dim pi As Double, D As Iteger, DF As Iteger Dim Y(5) As Double Dim a(, ) As Sigle, tau As Double Sub iit_wave() pi 4# * At(#) ' π 円周率 ' Frequecy (Hz) # / ' 周期 (sec) dt / 5 ' t t_max * 4 ' 積分範囲 DF 4 ' '------------------------------------------------------------- a.5 ' 初期値 a(, ).5: a(, ).: a(, 3) -.7: a(, 4) -. a(, ).8: a(, ) -.4: a(, 3).: a(, 4).3 '------------------------------------------------------------- t #: I Do I I + X(I) t Y(I) a For K o DF Y(I) Y(I) + a(, K) * Cos(K * # * pi * * t) + a(, K) * Si(K * # * pi * * t) Next K Cells(I, ) X(I) Cells(I, ) Y(I) t t + dt Loop While t < t_max + dt D I tau t_max / (D - ) Cells(, 3) D '**************** DF '**************** Ed Sub Sub DF() For N o 3 c #: s # N / D / tau For K o D - c c + Y(K) * tau * Cos(-# * pi * K * N / D) s s + Y(K) * tau * Si(-# * pi * K * N / D) Next K Cells(N, 4) Cells(N, 5) c Cells(N, 6) s Cells(N, 7) Sqr(c ^ + s ^ ) Next N Ed Sub
3.8 FF (FF: Fast Fourier rasorm) まず 最低何個の点で波を表現できるか考えてみる ω では最低 個のデータ店が必要になる ω 個 ω 4 個 ω3 6 個 ω4 8 個 ω5 個 図 3 データ点と波数 つまり データ点が N 個の場合 N/ 個の波数の波まで作れる これをサンプリング理論と呼ぶ サンプリ ング理論を考慮して周期のわからないデータを (3 4) 式の DF で計算すると 莫大な計算量になる 例 えば 周波数 4Hz までの音のデータをとる サンプリング理論から 秒間で N8 点のデータが必 要になる m 秒観測するならば 8 点のデータを取らなければならない (τ.5m 秒 ) ひとつの周波 数で 8 回 ( kτ) N π N τ e ik/ の和の計算を行う また 図 3 より 8 個のデータに対して 4 個 k の ω が作れる 結局 8 43, 回の計算をしないといけない そこで 計算回数を減らすために FF が考え出された DF の計算では もとの (t) に cos(ωt), cos(ωt), cos(3 ωt), cos(4ωt), cos(5ωt), をかけるが 図 3 を見 るとある時間 t で同じ値をかけていることがわかる この 同じ値の計算を省略すれば計算回数を減らせることができ πi/ N る (3 4) 式のτ 秒として W e を導入する と N N k G ( k) W (3 5) N k 図 3. cos(ωt) cos(3ωt) cos(5ωt) πi/ 8 つぎに N8 の場合の計算をすると W e となるので W 8 W, W 9 W の関係が求まる のとき G() ()W + ()W + ()W + (3)W + (4)W + (5)W + (6)W + (7)W のとき G(/8)()W + ()W + ()W + (3)W 3 + (4)W 4 + (5)W 5 + (6)W 6 + (7)W 7 のとき G(/8)()W + ()W + ()W 4 + (3)W 6 + (4)W 8 + (5)W + (6)W + (7)W 4 3 のとき G(3/8)()W + ()W 3 + ()W 6 + (3)W 9 + (4)W + (5)W 5 + (6)W 8 + (7)W (3.6) W 8 W, W 9 W の関係を利用すると 結局
() () () (3) (4) (5) (6) (7) G(/8) W W W W W W W W G(/8) W W W W 3 W 4 W 5 W 6 W 7 G(/8) W W W 4 W 6 W W W 4 W 6 G(3/8) W W 3 W 6 W W 4 W 7 W W 5 G(4/8) W W 4 W W 4 W W 4 W W 4 G(5/8) W W 5 W W 7 W 4 W W 6 W 3 G(6/8) W W 6 W 4 W W W 6 W 4 W G(7/8) W W 7 W 6 W 5 W 4 W 3 W W 表 3. これを偶数と奇数に分けると 偶数 () () (4) (6) 同じ G(/8) W W W W G(/8) W W W 4 W 6 G(/8) W W 4 W W 4 G(3/8) W W 6 W 4 W G(4/8) W W W W G(5/8) W W W 4 W 6 G(6/8) W W 4 W W 4 G(7/8) W W 6 W 4 W 奇数 表 3. () (3) (5) (7) G(/8) W W W W G(/8) W W 3 W 5 W 7 G(/8) W W 6 W W 6 G(3/8) W 3 W W 7 W 5 G(4/8) W 4 W 4 W 4 W 4 G(5/8) W 5 W 7 W W 3 G(6/8) W 6 W W 6 W G(7/8) W 7 W 5 W 3 W 奇数の場合は W をくくりだすと 表 3.3 のとき G(/8) odd W { ()W + (3)W + (5)W + (7)W } のとき G(/8) odd W { ()W + (3)W + (5)W 4 + (7)W 6 } のとき G(/8) odd W { ()W + (3)W 4 + (5)W + (7)W 4 } 3 のとき G(3/8) odd W 3 { ()W + (3)W 6 + (5)W 4 + (7)W } 4 のとき G(4/8) odd W 4 { ()W + (3)W + (5)W + (7)W } 5 のとき G(5/8) odd W { ()W + (3)W + (5)W 4 + (7)W 6 } 6 のとき G(6/8) odd W { ()W + (3)W 4 + (5)W + (7)W 4 } 7 のとき G(7/8) odd W 3 { ()W + (3)W 6 + (5)W 4 + (7)W }
W をくくりだした奇数は次のようにまとめられる () (3) (5) (7) 同じ G(/8) W W W W G(/8) W W W 4 W 6 G(/8) W W 4 W W 4 G(3/8) W W 6 W 4 W G(4/8) W W W W G(5/8) W W W 4 W 6 G(6/8) W W 4 W W 4 G(7/8) W W 6 W 4 W 表 3.4 偶数と奇数部分で上の部分と下の部分が共通なので計算回数が半分になることがわかる 但し偶数の部分で W を 4 回かけないといけないので 6 + (6+4) 36 回になる 64 回 36 回 同じ操作をするとさらに回数が減る (k) の偶数部分の (k) を p(k) とする k k となることに注意して ( W k V k ) W 8 W なので V 4 V となる p() p() p() p(3) G(/8) V V V V G(/8) V V V V 3 G(/8) V V V V G(3/8) V V 3 V V G(4/8) V V V V G(5/8) V V V V 3 G(6/8) V V V V G(7/8) V V 3 V V が求まる 偶数と奇数の部分に分けると 表 3.5 p() p() G(/8) V V G(/8) V V G(/8) V V G(3/8) V V G(4/8) V V G(5/8) V V G(6/8) V V G(7/8) V V 表 3.6 3
また 奇数部の V をくくりだしたものは V 4 V を考慮して p() p(3) G(/8) V V G(/8) V V G(/8) V V G(3/8) V V G(4/8) V V G(5/8) V V G(6/8) V V G(7/8) V V 表 3.7 問題 N8 の DF の計算 (64 回 ) は FF の計算で何回に減らすことができるか? FF の計算では 半分 半分 と分けて考えているので データは 個でなければならない 周期のわからない (t) を DF で計算したが ( 図 3 9) 同じ (t) を FF のプログラムで計算してみる G() 7 6 5 4 3 3 4 5 6 (Hz) 4
Dim pi As Double, D As Iteger, Poly As Iteger Dim Xr(5) As Double, Xi(5) As Double Dim ce(, ) As Sigle, tau As Double Sub iit_wave() pi 4# * At(#) ' π 円周率 ' Frequecy (Hz) # / ' 周期 (sec) tau / 5 ' τ (sec) D ^ ' umber o DAA '------------------------------------------------------------- Poly 4 a.5 ' 初期値 ce(, ).5: ce(, ).: ce(, 3) -.7: ce(, 4) -. ce(, ).8: ce(, ) -.4: ce(, 3).: ce(, 4).3 '------------------------------------------------------------- t # For I o D t (I - ) * tau Xr(I) a ' Real part (Iput DAA) Xi(I) # ' Imagiary part For K o Poly Xr(I) Xr(I) + ce(, K) * Cos(K * # * pi * * t) + ce(, K) * Si(K * # * pi * * t) Next K Cells(I, ) t Cells(I, ) Xr(I) t t + dt Next I '**************** FF_D '**************** For I o D / I / D / tau Cells(I, 4) Cells(I, 5) Xr(I) Cells(I, 6) Xi(I) Cells(I, 7) Sqr(Xr(I) ^ + Xi(I) ^ ) Next I Ed Sub ' 5
Sub FF_D() Dim s(5) As Double, c(5) As Double Dim A As Double, B As Double, da As Double Dim M As Iteger, H As Iteger, G As Iteger, P As Iteger, Q As Iteger M Log(D) / Log() ' N ^m A # da # * pi / D ' π/n For I o D / s(i) Si(A) ' W^k exp{ -i(π/n)k } c(i) Cos(A) A A + da Next I L D H For G o M L L / K For Q o H P For I K o L + K - J L + I A Xr(I) - Xr(J) B Xi(I) - Xi(J) Xr(I) Xr(I) + Xr(J) Xi(I) Xi(I) + Xi(J) I P he Xr(J) A Xi(J) B Else Xr(J) A * c(p) + B * s(p) Xi(J) B * c(p) - A * s(p) Ed I P P + H Next I K K + L + L Next Q H H + H Next G J D / For I o D - K D I J < I he Dummy Xr(I): Xr(I) Xr(J): Xr(J) Dummy Dummy Xi(I): Xi(I) Xi(J): Xi(J) Dummy Ed I Do K K / I J > K he J J - K Else Exit Do Loop J J + K Next I Ed Sub 6