管理人個人メモ/フーリエ変換周りの色々

Last-modified: 2010-09-27 (月) 17:05:20

参考文献

わかりやすいので追加
http://homepage2.nifty.com/takeuchiyosinori/CCP005.html

http://ufcpp.net/study/dsp/dft.html
が一番わかりやすい気がする。
http://laputa.cs.shinshu-u.ac.jp/~yizawa/InfSys1/basic/chap3/index.htm
は詳しいのだけど、ちょっと難しい。(分かった後の補足的な感じで)

http://www.aclab.esys.tsukuba.ac.jp/~ohbuchi/fourier/node1.html
は、わかりやすいところは神的にわかりやすいが、複素フーリエとかわかりにくいところもあり。

要点がよくまとまっているが、虫食いがある。
http://gijyutsu-keisan.com/mech/analysis/fourier_a/fourier_a_1.html#pagetop

忙しい人のためのフーリエ変換

1.http://www.kurims.kyoto-u.ac.jp/~ooura/fftman/ftmn1_2.html#sec1_2
  よりソースをコピペ。3つともFFTのソースだが、3つ目がオススメ
  ※HSPに移植する場合は、HSPでの#deffunc内の変数は
   C言語でいうところのstaticに相当し、再呼び出しの際に値が残っているので注意
   具体的に言うと、i = 0; とirev = 0; が必要だ。
2.ar[]に元データを入れる, nはデータの個数(2の累乗),
  thetaには、theta は ±2*PI/nを入れる(±どっちを入れても同じ?)
  ※ただし、HSPで2*M_PI/nとすると、0になってしまう。 2.0*M_PI/nとすべし
3.ar[]とai[]に結果が返ってくる
4.sqrt( ar[j]*ar[j] + ai[j]*ai[j] )がそのインデックスjにおける信号の大きさ
  ※間違えてar[]だけを見ないように注意する
5.インデックスjにおける周波数は、(分解能) x j である。
  分解能の計算式は...(執筆中)
6.インデックスjを0から始めてもよいが、低周波をはじきたい場合は、
  たとえば、15から始めれば、(分解能)x15以下の低周波を見ないようにすることができる。

フーリエ級数展開

直交

直交の関係にある関数の積を定積分すると値が0になる。
例:∫0~2π [sinx cosx] dx = 0
sinxとcosxは直交している

sin,cosの積(直交について)

周期の違うsinは直交
周期の違うcosも直交
sinとcosは、周期によらず直交

同じ周期の積の∫0~2πはπになる

y=acosx+bsinx

自然界では、0や1から始まるとは限らないので、
足し算して、グラフの位相を変化させる
(出発点の値を変えることができる)
※実際にグラフの形を考えて足し算してみよ。
※あるいは、半径aと半径bの円を考えて、acosxとbsinxを考えて、その合成を考えてみてもOK

位相をずらすのに、sin(x+θ)のようにしてはダメなのか?

計算が面倒くさくなる。

sinとcosだけで様々な位相を表せるといえる根拠は?

直交。
x軸、y軸の一次独立な組で任意の座標を表せるのと似ている。
難しく考えなくても、単位円状でacosθとbsinθを考えて、その合成(平行四辺形)を考えれば、
任意の位相を表せていることがわかる。
あるいは、三角関数の合成の公式を思い出してもOK

sinの合成の例

an sin nθ
でnを1~40、an=1,1/2,1/3,1/4,~1/40
として足し算すると鋸波になる。

また、同じ条件で、奇数だけ取り出して同じように足し算すると
矩形波になる。

フーリエ級数とは

F(x) = 1/2 a0 + Σn=1~∞ [an cos nx + bn sin nx]

フーリエ変換で一定の周期を持っていなければならない理由は?

フーリエ級数がsinとcosだけでできているから。
(フーリエ級数で作り出せるものが周期関数だけだから)
ちなみに、フーリエ級数で作られる周期は、一番、遅い周期のものに合致する。
(実際に波形を考えて足し算してみればわかる)

そもそもフーリエ変換ってなんだっけ?(今さら)

フーリエ級数の逆の考え方
フーリエ級数の係数anとbnがわかる⇒フーリエ級数の足し算の結果がわかる

↓逆

フーリエ級数の足し算の結果がわかる(原波形)⇒フーリエ級数のanとbnがわかる。
⇒このanとか、bnをフーリエ係数と呼ぶ

sin nxとか、cos nxのnが周波数に対応しているわけで、スペクトルがわかる


※厳密には、この手順はフーリエ級数展開とよぶ

やり方

1.波形を無理やり、ある時間にちょん切る
2.ちょん切った時間を(フーリエ級数の)最低周期とみなす
  ⇒たとえば、1秒区切りにしたら、1Hzってこと
3.「フィルター(後述)」を使って、a1,a2,...と1つずつフーリエ係数を抽出

フィルターって?

ヒント:直交=「異なる位相同士の」掛け算の定積分の結果が0
an cos nxだけを残す場合は、全体にcos nxを掛け算し0~2πで定積分
an cos^2 nx = □ ⇒ an = □/π
※∫0~2π cos^2 nx = π

つまり、
an = 1/π ∫0~2π F(x) cos nx dx
bn = 1/π ∫0~2π F(x) sin nx dx

残ったa0は?

ヒント:フーリエ級数がsinとcosであることに注目
    ⇒∫0~2πでみんな0になる
つまり、
a0 = 1/2π ∫0~2π F(x) dx

でも、実は、a0はcosのフーリエ級数で、a0だけ特別扱いしなくても、
前のanを求める式と同じになる。
ただし、値は2倍になってしまう⇒フーリエ級数の謎の1/2の理由はコレ。
要するに1/2をつけることでつじつまあわせをしているだけ。

で、どうやってスペクトルを得るの?

周波数がnに対応している波は
an cos nx + bn sin nx
で表されるから、
そこでの大きさはrn = √(an^2 + bn^2)
で、rnを左から順に並べればスペクトルの完成!

フーリエ級数展開とフーリエ変換

http://detail.chiebukuro.yahoo.co.jp/qa/question_detail/q1410197109
によると、
上記の手法はフーリエ級数展開に属するようだ。
==
フーリエ級数展開で得られるのは,ある領域内にある波形を基本周波数の整数倍の振動で表現したときの,各周波数成分の係数ですね。つまり,飛び飛びの離散変数としての係数を得るわけです
==
フーリエ変換は領域が限定されず,そのために基本周波数は無限小になるので変換後の係数ももはや離散的ではなく連続変数になります。計算自体も無限和ではなく積分となります。

複素フーリエ級数

フーリエ級数を指数関数で表す

参考http://www.akita-nct.jp/yamamoto/lecture/2006/3E/7th/html/node3.html
三角関数がイヤなので、無理やり指数関数を使って、フーリエ級数を表すもの。


オイラーの公式より、
cosx = e(ix) + e(-ix) / 2, sinx = e(ix) - e(-ix) / 2i
に注目して、フーリエ級数に突っ込む。

すると、
f(x) = a0/2 + Σn=1~∞[1/2(an-ibn)*e(inx) + 1/2(an+ibn)*e(-inx)]

※元のフーリエ級数からモノ自体は変えていない。
※また、右辺に複素数を含むが、キャンセルされてゼロになる。

=
ここでさらに、かなり無理があるが、以下のように定義する。
c0 = a0/2, cn = 1/2(an-ibn)*e(inx), c-n = 1/2(an+ibn)*e(-inx)
そうして、以下の複素フーリエ級数を得る
f(x) = Σn=-∞~∞[cn*e(inx)]

複素フーリエ級数の係数

c0 = a0/2 = (1/2π)∫-π~π[f(x)dx]
cn = 1/2(an-bn) = (1/2π)∫-π~π[cos nx - isin nx]dx

               = (1/2π)∫-π~π[e(-inx)]dx

c-n= 1/2(an+bn) = (1/2π)∫-π~π[cos nx + isin nx]dx

               = (1/2π)∫-π~π[e( inx)]dx

よくみると、3つまとめて、
cn = (1/2π)∫-π~π[e(-inx)]dx [n = 0, ±1, ±2,...]
とできる。

また、cn* = c-n

周期が2πでないとき

http://ufcpp.net/study/dsp/fourierseries.html
には、周期Tが2π以外の関数に関しては、変数tを2πt/T
で置き換えることにより、とある。ただの変数変換。
f(t)にt=2πt/Tを入れたから、f( 2πt/T )になるのではないか?f(t)のままなのは、
おかしいと思うかも知れないが、よく見ると、f(t)がt=2πt/Tが代入されたものに変わっている。
t=2πt/Tではなく、ωtで置き換えることもあるが、ω=2π/Tだから、結局同じ。
==
(リンク先が学内ページで、パーミッションを変えたのか?今はリンクできず)
http://www.tsunami.civil.tohoku.ac.jp/hokusai2/class/spec/02fourier.pdf
にもあった。これは
an = 1/π ∫0~2π F(x) cos nx dx
において、x = ωtとおいて、
an = ω/π∫0~(2π/ω)f(t)cosωt dt
となり、T=2π/ωであるから
an = 2/T∫(0~T)f(t)cosωt dt
⇒cosωtがcos(2πt/T)になっている場合もあり。
⇒0~Tの代わりに-T/2~T/2になっている場合もあり。
⇒また、周期がTであるから、an = 1/T∫(-T~T)f(t)cos(2πt/T)dtみたいになっている場合もあり。
てんでバラバラw

ちなみに、cos nx同士、sin mx同士の掛け算は、n≠mならば、0, n=mならば、π/ω
cosnxとsinmxの掛け算は0となり、直交は成立している。
変数変換しただけなのだから、当たり前か?
==
また、このとき、フーリエ級数は
f(t) = a0 + ∑n=1~∞( ancos2πnt/T + bnsin2πnt/T )
となる。
==
で、複素フーリエは、
cn = 1/T∫-T/2~T/2[ f(t)exp(-i2πnt/T )dt
f(t) = ∑n=-∞~∞[ cnexp( i2πnt / T )]
こいつも、∫区間が-T~Tになっていて、代わりに1/2Tになっていたりと色々バラバラw

フーリエ変換

フーリエ級数+複素フーリエ+周期が2πでないケース+拡張=フーリエ変換

http://laputa.cs.shinshu-u.ac.jp/~yizawa/InfSys1/basic/chap4/index.htm
http://www.geocities.jp/the_cloudy_heaven/laboratory/fouriertrans/fouriertrans.html
を参考
==
cn = 1/T∫-T/2~T/2[ f(t)exp(-inωt )]dt
に両辺Tをかけて
cn * T = ∫-T/2~T/2[ f(t)exp(-inωt )]dt
左辺をF(jnω)と置くと、
F(jnω)= ∫-T/2~T/2[ f(t)exp(-inωt )]dt
さらに、ω=2π/Tであるが、Tを無限大にすると、
ωは0に近づき、離散的なnωは連続量に変わる。
(なぜなら、n=-∞,...,-1,0,1,∞で、nω=-∞,...,-ω,0,ω,...,∞だけど、
 ωがクソ小さいと、nωの間隔もクソ小さくなるから。
 nωが離散値で、その離散間隔がωになっているのがミソ)


で結局、
F(jnω) = ∫-T/2~T/2[ f(t)exp(-inωt )]dt
で、
T⇒∞
nω(離散)⇒ω(連続)
として、
F(jω)=∫-∞~∞[ f(t)exp(-iωt)]dt
をフーリエ変換と呼ぶ。Tを無限大に飛ばすとは、周期が無限大ということだから、
フーリエ級数のときみたいに、擬似的に周期関数であるとみなさなくて良くなったわけだ。


逆変換は
f(t) = ∑n=-∞~∞[ cnexp( i2πnt / T )]
だったから、cn*T=F(jnω)⇒cn=F(jnω)/Tを考慮して
f(t) = ∑n=-∞~∞[ F(jnω)/T * exp( inωt ) ]
さらに、ω=2π/Tだから、×ω÷ω=ω*T/2πなので、
   = ∑n=-∞~∞[ ω÷T/2π*F(jnω)/T * exp( inωt ) ]
   = 1/2π * ∑n=-∞~∞[ ω * F(jnω) * exp( inωt ) ]
ここで、同様にTを無限大に飛ばすと、ω=2π/T ≒ 0 = dω(クソ小さいω)
nω(離散)⇒ω(連続となる)
そして、Σは連続量になり∫となる。
   = 1/2π * ∫-∞~∞[ dω * F(jω) * exp( iωt ) ]
   = 1/2π * ∫-∞~∞[ F(jω) * exp( iωt ) ]dω


※変換前の元信号がf(t)で、変換後がF(t)ね。

F(jω)=∫-∞~∞[ f(t)exp(-iωt)]dt
f(t)   = 1/2π * ∫-∞~∞[ F(jω) * exp( iωt ) ]dω

フーリエ変換とフーリエ級数の問題

http://nabe.blog.abk.nu/whats-fft
より
フーリエ変換は無限に続く連続信号に対するもので、フーリエ級数展開は有限長の連続信号に対するものです。どちらも、sin波とcos波の加算として表現されることから、その信号の持つ周波数ごとの信号の強さを知ることができます。

ですが、ここに問題があり、コンピューターが無限に続くものを扱えないのは当然としても、そもそも連続信号を扱うことができません。

(ディラックの)デルタ関数とは

※この項目は、離散フーリエを導く段階で語られることがあるが、
 離散フーリエの導出方法は複数あって、区分求積から導く流派ではデルタ関数が出てこないこともある。
原点以外は0だが、-∞~∞まで∫すると値が1になるというもの。
http://homepage2.nifty.com/eman/electromag/delta.htmlがわかりやすい
1x1の正方形(面積1)を横に縮めていくと、横が0になったときに、縦が無限大に飛ぶ。
連続性はない。(無限大を少しでもずれると0に飛ぶので明らかに連続でない)
こいつの存在意義は、リンク先が詳しい。
なぜ積分の結果が 1 になるようにしたかというと、参考サイトでは点電荷の積分の値を
正しくなるようにするために定義したもののようだ。
ようするに、関数の一部を取り出したいときに使えるということ。
δ(t-nT)だったら、t=nTのときだけ、(つまりδの中身が0のときだけ)1だけど、
それ以外は、0になるよという使い方ができる。

1と0を便利に扱うという点で、クロネッカーのデルタと似ているが、こいつは積分絡みで出てくるという理解でよさそう。
つまり、0の地点で∞なのはどうでもよくて、積分すると、その部分だけ1にできるということの方が重要なんだと思う。

窓関数

フーリエ変換を離散フーリエ変換に直したときに、一定の周期があることを考慮した。
だから、対象波形は周期があることが望ましい。
だが、任意の範囲を切り捨てると、周期がない状態になる。
切り捨てた部分の左端と右端が異なる=周期的でない⇒終端辺りに超高周波が含まれている扱いになり、フーリエ変換の結果にゴミが入る。ゴミがどの程度高周波になるかが読めないので、場合によってはデータをじゃまする可能性がある。
⇒ローパスをかけることが推奨されているが、それで切り取れるか否かは用途による。
⇒矩形で切り取るから、両端に急激なずれが生じるわけで、矩形で切り取らないで、滑らかに切り取るためにかけ合わせる関数を窓関数という。
※ちなみに、sinなどの周期関数であっても、切り取る場所が悪いと両端にずれが生じるのは同じこと。
http://gijyutsu-keisan.com/mech/analysis/fourier_a/fourier_a_1.html#pagetopによると
 この打ち切りによる誤差をリーケージエラーというらしい。言葉はどうでもいいが、検索には重要だよね。

窓関数の決定法

用途による。色々ある。面倒な場合は矩形でOK。

離散フーリエ変換とフーリエ級数の違い

http://www.kurims.kyoto-u.ac.jp/~ooura/fftman/ftmn1_1.html#sec1_1_3
離散フーリエ変換は通常の Fourier 変換の無限区間積分を有限の和で書き換えたもので,時間領域,周波数領域ともに離散化された Fourier 変換のことです
Fourier 級数展開は,周波数領域でのみ離散化された変換に相当します

その他

言葉を認識するのは難しいらしい。特に、子音が音が短くて取りにくいらしい。
音程は結構正確にいける。リコーダなどはsin波に近くて、とりやすい。
認識させる側(人間)もうまく認識させられるように学習できるらしい。