フーリエ級数展開をExcelのFFTでシミュレートする

たとえば、連続関数 f(x) = x(π-x) を３次まで展開すると

　-(π/3)(-6sin(x)+3sin(2x)-2sin(3x)+π)+4cos(x)-cos(2x)+(4/9)cos(3x) ……(#)

になるようです。wolframa で確認した限りでは、３次ではあまりいい近似になりません。wolframa のフリー版では4次以上も計算してくれるようですが、グラフまでは示してくれません。

　で、この
　　「あまりいい近似でない３次のフーリエ級数展開」
を、Excel の FFT を使って再現するにはどうすればいいのかというのが質問です。

　まずf(x)から128個切り出した
　　g[n](0≦n≦127)
を用意する必要がありますが、元の関数f(x)の定義域にだいたい合わせて
　　g[t] = g[t(n)](0 ≦ t(n) ≦ 3.175 : 0≦n≦127 : t = t(n) の刻み幅は 0.05)
　　t(0) = 0, t(1) = 0.05, t(127) = 3.175
とします。このg[t]に対しFFTを実行し、続けてIFFTを実行すれば当然元のg[t]が回復されます。これは元の関数でいえば f(x) = x(π-x) を128次まで展開することに相当すると思われます（違うのだろうか？）。いま、FFT で実現したいのは f(x) = x(π-x) を３次まで展開することに相当することです。

　g[t]のFFTの結果からf(x)の３次に相当すると思われる

　　　g[t(1)], g[t(2)], g[t(3)]
　　　g[t(125)], g[t(126)], g[t(127)]

と直流成分 g[t(0)] 以外の値をすべて削除して（つまり 0 にして）IFFT を実行したのですが、これでは全然だめです。f(x) = x(π-x) の３次よりずっといい近似になってしまうのです（笑）。何か決定的な勘違いをしていると思うので、ご教示いただけたら幸いです。
参考画像
https://imepic.jp/20230703/729270

※ExcelのFFTは使いづらいので、ネットで拾ったVBAコードを利用しています。

No.5 ベストアンサー 回答者： stomachman 回答日時： 2023/07/04 04:25

No.1へのコメントについて。

> 反応しませんでした。

FFTのcommandがどうなってるんだかはdocumentをガッツリ読まないと。

　ともあれ、手直しの仕方について。もしFFTを使ってx(π-x) を0～πの区間でフーリエ級数展開したいのなら、

> これは元の関数でいえば …（違うのだろうか？）。

　違うっす。刻み幅0.05、なんてやりかたではうまくいきません。
●第一に、x(π-x)というのを
　　f(x) = z(π-z), where z = x - [x/π]π （[t]はtを越えない最大の整数（gaussの記号））
という周期πの関数として扱う必要がある。
●第二に、N点のFFTは周期Nの周期関数だけを相手にしますので、周期がNになるように変換する必要がある。
　だからこの場合には
　　g(n) = (πn/N)(π - (πn/N))　　(n=0,1,....,N-1)
を使わなくては。

　一方、wolframaの方もこれではダメみたいですね。（皆さんの回答を拝見すると、少なくともデフォルトでは） -π≦x＜πを周期とする関数を展開するみたいですから、変数変換して周期を2πにする。さらにx≦0のところが正しく周期通りになるようにする。FFTに与える関数のグラフで言えば、「n=N/2〜N-1の部分を切り取って、それをn<0の部分に平行移動した形」を作るわけで、
　　h(x) = |x|(2π - |x|)
のフーリエ級数展開をやらせなくちゃいけません。-π≦x＜πでプロットすれば x=0を中心として左右対称の「Υ」のような形になるわけです。

　ところで、こういう簡単な関数なら、じかに積分して確認できるでしょう。
　　∫[-π〜π] h(x) e^(-2inx) dx
　　= 2∫[0～π] h(x) cos(2nx) dx
　　= 2∫[0～π] x(2π-x) cos(2nx) dx
　　= if n=0 then 4(π^3)/3, else -π/(n^2)
というわけで、展開の係数は実数、つまりcosine成分だけです。

　なお、両者はそれぞれの変数変換に伴うヤコビアン（定数の係数）を（手抜きして）無視したせいで、互いに定数倍の違いが生じます。が、ま、ややっこしいことを考えるより、n=0での値が一致するように規格化した上で比較すれば良いですね。

この回答へのお礼

丁寧な回答まことにありがとうございました。深く感謝いたします。

お礼日時：2023/07/04 10:42

No.4 回答者： eatern27 回答日時： 2023/07/04 00:19

https://reference.wolfram.com/language/ref/Fouri …

を見る限り、デフォルトでは積分範囲は-πからπで、
積分区間の幅を2πから変えるにはFourierParametersを設定すれば良いようです。

積分区間の中心を0からずらすパラメータに関する記載はないので、t→t-π/2と変数変換して、積分区間をπ/2からπ/2にした式を入力するしかないかもしれませんね。wolframの関数はほぼ知らないので、他の関数で計算できるかもしれませんがそこまでは調べてません。

No.3 回答者： Tacosan 回答日時： 2023/07/03 23:50

例から推測するに, wolframalpha では

-π～π
でフーリエ級数を計算してるんだと思う. 関数 x(π-x) をこの範囲で周期化すると x= ±π で連続しないので, 結果境界ふきんで値があばれる. Gibbs 現象の一種だね.

境界で値がつながるように
(π+x)(π-x)
とすればそこそこ見れた形になる. 境界のところは (導関数が連続じゃないので) やっぱりくるしいけど.

No.2 回答者： eatern27 回答日時： 2023/07/03 23:38

>　-(π/3)(-6sin(x)+3sin(2x)-2sin(3x)+π)+4cos(x)-cos(2x)+(4/9)cos(3x) ……(#)

これは具体的に何を計算して出てきたのですか？

それ次第だとは思うけど、何かしらのツールを使って計算していて、デフォルトのパラメータを使っている部分があるのなら、
>0 ≦ t(n) ≦ 3.175
これに対応する部分を0 ≦ t(n) ≦ 2πまたは-π≦ t(n) ≦ πだと思えば、良いかも？

この回答へのお礼

回答まことにありがとうございます。
>　-(π/3)(-6sin(x)+3sin(2x)-2sin(3x)+π)+4cos(x)-cos(2x)+(4/9)cos(3x) ……(#)
>　これは具体的に何を計算して出てきたのですか？

　wolframa で

FourierSeries[x(π-x),x,3]

と入力しました。定義域を 0≦x≦π と指定したいのですが、できません。

FourierSeries[x(π-x),x,3], {0,π}

でもダメでした。

お礼日時：2023/07/03 23:46

No.1 回答者： stomachman 回答日時： 2023/07/03 22:23

FFTがやってるのは、x=0〜127の範囲の波形を無限に繰り返した周期関数

　　f(x) = f(x+127n) （ただしnは任意の整数）
の離散フーリエ級数展開です。ご質問の場合ですと、f(x) = f(128-x)という格好をしているから、sine成分は全て0で、cosineの項だけが残るはず。（複素FFTを使ったのなら、すべての係数が実数になる。）
　wolframaではおそらく、「x=0〜127の範囲以外では上記の周期関数とは似ても似つかない、何か別の関数」のフーリエ級数展開をやったんでしょう。

この回答へのお礼

回答まことにありがとうございす。

> x=0～127の範囲以外では上記の周期関数とは似ても似つかない、何か
> 別の関数のフーリエ級数展開をやったんでしょう。
　ああ、なるほど。であれば、wolframaでフーリエ級数展開するとき、定義域を指定するにはどうしたらいいのでしょうか？

FourierSeries[x(π-x),x,3]

では図のような結果になります。

FourierSeries[x(π-x),x,3],(0≦x≦π)

としたら、反応しませんでした。

お礼日時：2023/07/03 23:16