関数のフーリエ展開の数値的な扱い

質問

f(x)は |x| < a では正の有限値を持つ滑らかな関数で、
|x| > a ではゼロとします。x = a は連続につながっています。
たとえば
　f(x) = cos(x) (|x|≦π/2) , 0 (|x|> π/2)
のような関数です。

この関数を間隔 L > 0 で並べた関数をφ(x)とします。
　φ(x) = Σ[n] f(x + nL)

L が a に比べて十分に大きいような場合、すなわち
有限値を持つ山が、値がゼロである領域の谷をはさんで
飛び飛びに並んでいるような場合について考えます。

適当な基底関数{Vn(x)}で展開して
　φ(x) = Σ[n] Cn Vn(x)
という形にして、数値的に扱いたいのですが
三角関数系{exp(inx)}を基底にすると、かなり大きな n まで
和をとらないと精度が保てません。このような場合に、
少ない項数で精度よく展開できるような方法は
有りませんでしょうか？アイデアを頂けるとありがたいです。
よろしくお願いします。

回答No.2

stomachman · ベストアンサー

「精度が保てない」と仰る意味がはっきりしませんが…

[1] ご質問のような周期2πの周期関数f(x)の半周期分だけを取り出したものをg(x)とすると、矩形関数
B(x) = |x|＜π/2のとき1、さもなくば0
を使って
g(x) = f(x)B(x)
とすれば良い。逆にg(x)からf(x)を作るには（Σをn=0, ±1, ±2, ....　の総和として）
f(x) = Σ(g(x-2nπ))
= g(x)＊Σ(δ(x-2nπ))
とすれば良い。（ここで＊は畳み込み(convolution)、δはディラックのデルタ関数(本当は超関数)です。）

　矩形関数Bをかけ算することで切り取られる縁の持つ高周波成分がいろいろ問題になるという場合には、縁を持っている基底を使うのも一つの手です。たとえばウォルシュ・アダマール変換とか、wavelet変換の中でとんがった基底を持つ物を使うとか。しかし、その基底がf(x)とよほど相性が良い（うんと少ない項数で表せる）というのでない限り、計算量の観点から言うと、普通の（離散）フーリエ変換の方が実用的です。なにしろFFTという優れたアルゴリズムがあるんで、項数を多くしてもうんと速く計算できるから。

[2] というわけで、（お分かりなのかもしれないけれど、）普通のフーリエ変換の話。
　めんどくさい定数倍とか、横軸の単位とか、細かい話は抜きでやります。（フーリエ変換の定義の仕方（流儀）によって細かい違いが出ますので。）なので、以下はご自分で計算しなおしてくださいね。
　f(x) の三角関数による（普通の）フーリエ変換F(ω)は、コンボルーション定理により、
F(ω) = G(ω) (Σ(δ(x-2nπ))のフーリエ変換）
であり、
(Σ(δ(x-2nπ))のフーリエ変換）= Σ(δ(ω-m))
となりますから、結局
F(ω) = Σ(G(ω)δ(ω-m))
となる。δ関数に付く係数をFmと書くことにすれば、
F(ω)= Σ(Fm δ(ω-m))
Fm = G(m)
で、この係数Fmこそがフーリエ級数に他なりません。（つまり、フーリエ級数とは周期関数のフーリエ変換とほとんど同じこと。）
　さて、ご質問で例として挙がっているcosine関数の半周期分を取り出すには、
g(x) = B(x)cos(ax)
とすればいいですね。その三角関数による（普通の）フーリエ変換G(ω)は
G(ω) = (B(ax)のフーリエ変換)＊(cos(ax)のフーリエ変換)
となります。で、
(B(ax)のフーリエ変換) = sinc(ω/a) (ただし、sinc(x) = (x≠0のときsin(x)/x 、x=0のとき1）)
(cos(ax)のフーリエ変換) = δ(1/a)+δ(-1/a)
となるから、
G(ω) = sinc((ω+1)/a)+sinc((ω-1)/a)
である。
　もうちょっと一般的に、ある関数h(x)があって、
g(x) = B(x)h(x)
と書け、h(x)は滑らかで|x|→∞で指数関数的に0に漸近する場合を考えます。h(x)のフーリエ変換をH(ω)として
G(ω) = sinc(ω/a) ＊H(ω)
です。で、最終的に欲しいのは
Fm = G(m)
です。
　このFをどう使うのか、応用によってはHだけ分かっていればうまく計算が進められる場合だってあります。

お礼

丁寧なご回答ありがとうございました。

> ウォルシュ・アダマール変換とか、wavelet変換

これらについて調べてみましたが、ご指摘のとおり、
基底を変えて多少メリットがあったとしても、
FFTのメリットをしのぐ程では無いようですね。
勉強になりました。

回答No.3

stomachman

ANo.2のstomachmanですが、ANo.1に付けられたコメントを拝見しました。

　えとですね、段差のある波形をフーリエ級数で表しておいて、そのフーリエ級数を有限項で打ち切った和を作ると、元の波形がほぼ再現されるのだけれども、段差の縁に細いツノが生えます。どんなに項数を多くしてもそれが有限個である限り、ツノは（細くなるけれども）結構大きな高さを保って、決して小さくならない。この事には「ギブスの現象」という名前がついてます。
　こいつを消すには、段差を滑らかな斜面で置き換えるために、元の波形に平滑化を掛けておく（たとえばガウス関数を畳み込む）必要があります。（もちろん、Ckに重みを付けてやっても同じ事です。）

回答No.1

spring135

目的がよくわからないので見当はずれかもしれませんが、とにかく周期関数をスペクトルに分解したいならばＦＦＴが文句なしに答えを出してくれます。
http://ja.wikipedia.org/wiki/%E9%AB%98%E9%80%9F%E3%83%95%E3%83%BC%E3%83%AA%E3%82%A8%E5%A4%89%E6%8F%9B
またウェーブレット変換という方法もありますが使ったことはありません。

補足

早速のご回答ありがとうございます。
質問の意図が自分で読んでもわかり辛いので補足をさせて下さい。

目的は、大仰に言えば固体の電子のシミュレーションです。
電子の分布関数φ(r) = Σ[k] Ck exp(ikr)と仮定して
シュレディンガー方程式に代入、数値的にCkを求め、
求まったφを利用して物理量を概算します。

理想的な無限結晶(1)なら、滑らかな周期系なので
展開係数Ckは急速にゼロに収束します。Σの足し算を
kが小さいうちに打ち切っても、問題が生じません。
しかし、固体の中に大きな空洞があるような場合(2)、
Σの足し算を早めに打ち切ると、まともな物理量が
出せなくなってしまいます（電子密度がマイナスになる等）。
項数を多くしていくと、漸近的にまともな値に収束していきます。

(1)
… Na Cl Na Cl …
… Cl Na Cl Na …
(2)
… Na Cl　　　　　　　　　　 　　　Na Cl …
… Cl Na 　　　　　　　　　　　　　Cl Na …

素朴考えですが、こういう状態にぴったりの基底があればなぁ
と思ったので質問させていただきました。

でも、実際に計算する事を考えると普通にFFTがベストのようですね。
ウェーブレット変換についてはまた調べてみます。