周波数を抜きだそう

[size=150][b]このワークシートは[url=https://www.geogebra.org/m/twxxx3yq]Math by Code[/url]の一部です。[br][br]アプレット、背景、実装の順に見ていきましょう。[br][/b][/size]
フーリエ級数と周波数表示
時間領域関数f(t)と周波数領域関数F(ω)の連動
フーリエ級数展開は[br]特定区間の関数や周期関数を倍音周波数の純音に分解することだった。[br][br]今回は、[b]一般の関数を周波数を抜き出すフーリエ変換[/b]を[br]みていこう。[br][br]フーリエ変換がいろんな音に使えるとすごく便利になる。[br]たとえば、ギターやマイクをPCにつないで、オーディオデータとしてPC入力したとしよう。[br]その[color=#0000ff]データをフーリエ変換すれば、音の周波数、つまり音程がわかるはず[/color]だね。[br]つまり、その音程をMIDIデータとして活用できるよいうことだ。[br][br]最近のDAW(digital audio workstation)という[b]DTMソフト[/b]ならば、オーディオデータを[br]midiデータに変換してくれる機能がある。これで、midiデータをポチポチ入力しなくても、[br]すでにある音楽データをmidiデータにしてくれるということになる。[br][br]この機能はまさに、[b]フーリエ変換[/b]が基本にあるのではないでしょうか。
1.背景
フーリエ級数展開を周波数分析という視点で振り返ることからはじめよう。[br]・純音はy=a sin[ 2π f t] (f=440でA4の音)だった。[br]・ノコギリ波は[ーπ、+π]の区間で、f(t)=t, T=2πのときのフーリエ級数展開は[br]f(t) = 2 Σ(-1)[sup](n+1)[/sup]/n [sub][/sub] sin(nt) [math]=2(sint-\frac{1}{2}sin2t+\frac{1}{3}sin3t-\frac{1}{4}sin4t+...........)[/math] [br]wa(x)=Σ[(1/m) sin(m 220 frq*2 π x) (-1)^(m+1)]という形で実装した。[br]・矩形波は[0, π]区間で作れたら奇関数で、f(t)=1, T=2πのときのフーリエ級数展開は[br]f(t) = Σ(1-(-1)[sup]n[/sup]) /nπ [sub][/sub] sin(nt) [math]=\frac{2}{\pi}(sint+\frac{1}{3}sin3t+\frac{1}{5}sin5t+..........)[/math][br]wa(x)=Σ[(1/(2p+1)) sin((2p+1) 220 frq*2 π x)]という形で実装した。[br]・三角波はf(t) = t (0からπ/2) ;π-t (π/2, π)を偶関数化したフーリエ級数展開は[br]f(t)=π/4- 2/π{ cos 2t + 1/3[sup]2[/sup] cos 6t +1/ 5[sup]2[/sup] cos 10t + ......1/ (2m+1)[sup]2[/sup] cos2(2m+1)+.......}[br][math]=\frac{\pi}{4}-\frac{2}{\pi}\left\{cos2t+\frac{1}{3^2}cos3\cdot2t+\frac{1}{5^2}cos5\cdot2t+........+\frac{1}{\left(2m+1\right)^2}cos\left(2m+1\right)2t+......\right\}[/math][br]wa(x)=(π/4)+Σ[(-2/π)*1/((2p+1)^2)) cos((2p+1) 220 frq *2 π x ]という形で実装した。[br]これらをまとめてみよう。[br]フーリエ級数展開では、周期関数f(t)と倍音周波数の純音に和分解していることになる。[br]たとえば、基本周波数をfとすると、(周波数の倍率,振幅の倍率)のペアのリストができる。[br]・純音   ;(1,1)[br]・ノコギリ波;(1,1),(2,-1/2), (3, 1/3),(4,-1/4),(5,1/5),......................, ( k, (-1)[sup]k+1[/sup]/k ),.....[br]・矩形波  ;(1,1)    , (3, 1/3) ,(5,1/5),......................., ( k, (1-(-1)[sup]k[/sup])/k ),.....[br]・三角波  ;(1,1),    ,(3, 1/3[sup]2[/sup]) ,(5,1/5[sup]2[/sup]),,..................., (2k-1, 1/(2k-1)[sup]2[/sup]),.....[br][b][color=#0000ff][size=150]この数列や数式を最初のアプレットで実感してみよう。[br][/size][/color][/b][br]これを作ることができたのはフーリエ級数展開の公式のおかげだったね。[br][br][b][size=150]<フーリエ級数から複素フーリエ級数へ>[br][/size][/b]一般に、周期Tのf(t)のフーリエ級数はcos波やsin波のような周期関数の無限個のたし算だったね。[br][size=150]f(t)=a0+ Σ( a[size=150][sub]n[/sub][/size] cos[ 2π n t/T ] + b[size=150][sub]n[/sub][/size] sin[ 2π n t/T] ) [br]a0 = 1/T ∫ f(t) dt (-T/2, T/2) [br]a[size=150][sub]n[/sub][/size] = 2/T ∫ f(t) cos[ 2π n t/T ] dt ( -T/2, T/2) 、[br]b[size=150][sub]n[/sub][/size] = 2/T ∫ f(t) sin[ 2π n t/T ] dt (-T/2, T/2) [br][/size]のように。[br][br]・ここで、オイラー公式e[sup]it[/sup]=cos t+ i sin t, e[sup]-it[/sup]=cos t - i sin tを思いだそう。[br] e[sup]ix[/sup]はcosx, sinxで作れれているのだから、周期関数と言える。[br] cos t =( e[sup]it[/sup] +e[sup]-it[/sup])/2 , sin t =( e[sup]it[/sup] -e[sup]-it[/sup])/2i におきかえられるね。[br] さらに、k[sub]n[/sub]=2πn/T, x = k[sub]n[/sub] t , [b]c[sub]n[/sub]=(an-ibn)/2 (cn*=(an+ibn)/2)とすれば、[/b][br][b]sinx,cosxではなく、e[sup]ix[/sup]を基本の波として無限個たし算する複素フーリエ級数[/b]公式ができるはずだ。[br]knは単位円での回転で2πを1周としたときの周期をTとするとき、k1=2π/Tの何倍の[b][color=#0000ff]角周波数[/color][/b]かを表す。[br][br][size=150][b]f(t)[/b]=a0+ Σ( a[size=150][sub]n[/sub][/size] cos knt + b[size=150][sub]n[/sub][/size] sin knt )[br][size=150]=a0+1/2Σ [ a[size=150][sub]n[/sub][/size]( e[sup]iknt[/sup] +e[sup]-iknt[/sup])- ib[size=150][sub]n[/sub][/size](e[sup]iknt[/sup] -e[sup]-iknt[/sup])](1,∞)[br][/size]  =a0+1/2Σ [ (a[sub]n[/sub]-ib[sub]n[/sub])e[sup]iknt[/sup] +(a[sub]n[/sub]+ib[sub]n[/sub])e[sup]-iknt[/sup]](1,∞) [br]  =a0+Σ [ cn e[sup]iknt[/sup] +cn*e[sup]-iknt[/sup]](1,∞) [br]これから、n=0のときのa0も含むものとして、[br][b][color=#0000ff]複素フーリエ級数f(t) [br]=...+c[sub]-1[/sub]e[sup]-ik-1t[/sup]+c[sub]0[/sub]+c[sub]1[/sub]e[sup]ik1t[/sup]+......[br]= Σc[sub]n[/sub] e[sup]iknt[/sup](‐∞,∞)[br][/color][/b][size=100]cn=(an-ibn)/2をまとめると、[br]cn={an-i bn}/2[br] =[ 2/T ∫ f(t) cos[ 2π n t/T ] dt ( -T/2, T/2) -i 2/T ∫ f(t) sin[ 2π n t/T ] dt (-T/2, T/2) }/2 [br]cn= 2/T ∫ f(t) 1/2(cosknt- i sinknt) dt (-T/2, T/2) となるので、[br][/size][b][color=#0000ff][sub]複素フーリエ係数cn[/sub]=1/T ∫ f(t) e[/color][sup]-iknt[/sup][color=#0000ff] dt( -T/2, T/2)[br][/color][/b][/size]sin,cosにばらけていたものが、e[sup]ix[/sup]に一本化されてスッキリしたね。[br][b][color=#0000ff]周期Tの波の記述をカンタンにした、[/color][/b]ともいえる。[br][br][b][size=150]<複素フーリエ級数からフーリエ変換へ>[br][/size][/b]さらに、複素フーリエ級数の周期T→∞にしたものがフーリエ変換と逆フーリエ変換だ。[br]区間(周期)Tの関数f(t)[b][color=#0000ff] =Σc[sub]n[/sub] e[sup]iknt[/sup](‐∞,∞)のT→∞にすると[br]角周波数kn=2πn/T→0は連続化してkとなる。[br][/color][/b]複素フーリエ係数cnのT倍をFT(k)= ∫ f(t) e[sup]-ikt[/sup] dt( -T/2, T/2)とかこう。級数f(t)のうち、角周波数k[sub]n[/sub]の項はc[sub]n[/sub]e[sup]iknt[/sup]=1/T FT(kn)e[sup]iknt[/sup] = 1/2π *(2π/T FT(kn)e[sup]iknt[/sup]) とかけるね。[br]( )の中は幅2π/T=k1,高さFT(kn)eikntの長方形の面積で、その和がf(t)とみることができる。[br]だから、T→∞のとき、ωn→0で連続化する。すると、積分の意味から、[br]f(t)=1/2π(長方形の細分の和)=1/2π(dk * FT(kn)e[sup]iknt[/sup]の和)=[b]1/2π [/b][b]∫F(k) e [sup]ikt[/sup] dk (‐∞,∞)[/b] となる。[br]ここで、F(k) はT→∞だから、FT(k)=∫ f(t) e[sup]-ikt[/sup] dt[b]( -T/2, T/2) [/b]→ F(K)=∫ f(t) e[sup]-ikt[/sup] dt[b]( -∞, ∞) で求められる。[/b][br][b][color=#0000ff][size=150]つまり、f(t)のフーリエ変換はF(k)=∫ f(t) e[sup] -ikt[/sup] dt (‐∞,∞) [br]F(k)の逆フーリエ変換f(t)= 1/2π∫F(k) e [sup]ikt[/sup] dk (‐∞,∞) [br][/size][/color][/b][br]f(t)という時間領域の非周期関数から、F(k)という周波数領域の連続関数への変換[br]これが、フーリエ変換だということがわかったね。[br][color=#0000ff](例)[br][/color]時間領域の関数f(t)=1 ([math]-\frac{\varepsilon}{2}\le t\le\frac{\varepsilon}{2}[/math] )に対するフーリエ変換は[br] F(x)=∫f(x)e[sup]-ixt [/sup]dt(-∞,∞)=[e[sup]-ixt[/sup]/-ix](-ε/2, ε/2)=-1/ix(e[sup]-ixε/2[/sup]-e[sup]ixε/2[/sup]) [br]=-1/ix[ cos(xε/2)-i sin(xε/2)- (cos(xε/2)+i sin(xε/2)) ] = 2i sin(x ε/2)/ix = ε sin(x ε/2)/(x ε/2)[br]εを大きくすると、時間領域tが広がる。それに連動して、周波数はx=0に集中します。[br]εを小さくすると、時間領域tが狭くなり、連動して、周波数領域は広がるという逆の関係が[br]見られるところが面白いね。[br]これは、冒頭の2番目のアプレットで確認しよう。[br]
2.実装
[color=#9900ff][b][u][size=150]質問:フーリエ変換のイメージをgogebraで視覚化するにはどうしたらよいでしょうか。[br][/size][/u][/b][/color][br]フーリエ変換は、f(t)という時間領域の非周期関数から、F(k)という周波数領域の連続関数への変換[br]だった。[br]変換のイメージをつかむだけなら、f(t)が周期関数で、F(k)が離散値でもよいね。[br]それが、冒頭のアプレットだ。[br][br]このアプレットでは、[br][br]たとえばノコギリ波の時間関数は、[br]para=Sequence(n)という自然数列に連動して、[br]sw(x)=Sum(Zip(((1)/(m)) sin(m f*2 π x) (-1)^(m+1),m,para))[br]これをグラフィックビューで表示する。[br][br]周波数の分布リストは、[br]saw=Sequence($point(k,(((-1)^(k+1))/(k))),k,1,n)[br]別のグラフィックビューで表示する。[br]これを上下に並べてみると、連動がよくかるね。[br][br]他の三角波、矩形波でもどうようにできる。[br]それをスライダーをつけて、整数値を変化させることで、表示の切り替えができるね。[br]周波数の深さnを変えると、関数のSumの深さも点列の数も連動して変わえられるね。[br]

Information: 周波数を抜きだそう