フーリェ変換（特に一次元のFFT）

質問

フーリェ変換について勉強しているものです。このページ（http://www.kurims.kyoto-u.ac.jp/~ooura/fftman/ftmndl.html）で配布している、fft.tgz (71 KB)というファイルのdfst()という関数を使おうと思っているんですが、その関数の返り値（返り配列と言った正しい？）の詳しい意味がよくわからないのです。
---------------------------------------------------------
以下の変換は :
S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n
このように分割される :
S[2*k] = sum_j=1^n/2-1 (a[j]-a[n-j])*sin(pi*j*k/(n/2)),
◎ S[2*k+1] = sum'_j=1^n/2 (a[j]+a[n-j])*sin(pi*j*(k+1/2)/(n/2))

注意:sum_■^▲はΣを表しており、■から▲までの部分和となっています。
nは要素数です。
------------------------------------------------------
この◎の関数はどのような意味を持つのでしょうか？

回答No.1

stomachman · ベストアンサー

実務でお困りなのだろうと推測します。

[1]dfst()なる関数がどこにあるんだか、また、ご質問の式がどこから出て来たものだか分かりませんけれども、周期2nを持つ実奇関数a[j]の離散Fourier変換（しかも、nが偶数の場合）の計算方法の話ですね。

> この◎の関数はどのような意味を持つのでしょうか？

に直裁に回答するなら、「関数aの離散フーリエ変換の、sine成分のうち、2k+1番目の成分の振幅」ですが、ただしご質問の◎の式は(a[n/2]=0である場合を除いて)間違ってます。

[2]a[j]が周期2nを持つということは
a[j+2*n]=a[j]
という意味であり、a[j]が奇関数であるとは
a[-j]=-a[j]
を満たすということです。従って、さらに
a[j]=-a[2*n-j]
a[n-j]=-a[n+j]
も成り立つことが分かります。

[3]周期2nを持つ複素関数a[j]の離散Fourier変換をF[k]とすると、0≦k＜2nについて
F[k] = c[k] + i*s[k] (iは虚数単位)
c[k] = sum_j=0^2*n-1 a[j]*cos(2*pi*j*k/(2*n))
s[k] = sum_j=0^2*n-1 a[j]*sin(2*pi*j*k/(2*n))
です。
　これをa[j]が奇関数である場合に限定すると、計算するまでもなく
c[k]=0, 0≦k＜2*n
である。
　一方、a[j]が実関数である場合に限定すると、
c[2*n-k]=c[k]
s[2*n-k]=-s[k]
である。
　従って、a[j]が実奇関数である場合に限定すると、
c[k]=0
s[2*n-k]=-s[k]
である。さらに
s[0]=0
である。だから0＜k＜nについてs[k]を計算するだけでF[k]が分かることになります。

[4]そこでs[k]の計算を、a[j]の性質（s[2*n-k]=-s[k]）を使って簡単にします。
s[k] = sum_j=0^2*n-1 a[j]*sin(pi*j*k/n)
において、j=0のときとj=nのときにはsin(pi*j*k/n)=0なので、
s[k] = {sum_j=1^n-1 a[j]*sin(pi*j*k/n)}+{sum_j=n+1^2*n-1 a[j]*sin(pi*j*k/n)}
右辺の第二項を s[2*n-k]=-s[k] を使って書き換えると、
第二項 = sum_j=1^n-1 a[j]*sin(pi*j*k/n)
ところがこれは第一項とピッタリ同じなので、結局
s[k] = 2*{sum_j=1^n-1 a[j]*sin(pi*j*k/n)}
です。だから、ご質問の最初の式
S[k]=sum_j=1^n-1 a[j]*sin(pi*j*k/n)
は、s[k]の1/2倍である。

[5]ようやく本題です。nが偶数であるという条件を使って、S[k]の計算をさらに簡単にしようと試みます。
[5-1] kが偶数(k=2m)の場合、
j=n/2のとき　sin(2*pi*(n/2)*m/n)=sin(pi*m)= 0
だから
S[2m]={sum_j=1^n/2-1 a[j]*sin(2*pi*j*m/n)}+{sum_j=n/2+1^n-1 a[j]*sin(2*pi*j*m/n)}
です。右辺の第二項は
-sum_j=1^n/2-1 a[n-j]*sin(2*pi*j*m/n)
と書き換える事ができ、従って
S[2*m] = sum_j=1^n/2-1 (a[j]-a[n-j])*sin(2*pi*j*m/n)

[5-2] kが奇数(k=2m+1)の場合、
S[2*m+1] = sum_j=1^n-1 a[j]*sin(pi*j*(2*m+1)/n)
= {sum_j=1^n/2-1 a[j]*sin(pi*j*(2*m+1)/n)}+a[n/2]*sin(pi*(n/2)*(2*m+1)/n)+{sum_j=n/2+1^n-1 a[j]*sin(pi*j*(2*m+1)/n)}
右辺の第三項は
sum_j=1^n/2-1 a[n-j]*sin(pi*j*(2*m+1)/n)
と書き換えることができ、また第二項は
a[n/2]*sin(pi*(n/2)*(2*m+1)/n)=a[n/2]*sin((m+1/2)*pi)= mが偶数→a[n/2]、mが奇数→-a[n/2]
です。以上から、
S[2*m+1] = {sum_j=1^n/2-1 (a[j]+a[n-j])*sin(pi*j*(2*m+1)/n)}+a[n/2]*sin((m+1/2)*pi)
あるいは、
S[2*m+1] = {sum_j=1^n/2 (a[j]+a[n-j])*sin(pi*j*(2*m+1)/n)}-a[n/2]*sin((m+1/2)*pi)
と書いても同じです。
　だから、a[n/2]≠0のとき、◎の式は間違っています。

[6] この間違いを修正するには、
p[j]=a[j]-a[n-j] (0＜j＜n/2)
q[j]=a[j]+a[n-j] (0＜j＜n/2)
q[n/2]=a[n/2]
として、
S[2*m] = sum_j=1^n/2-1 p[j]*sin(2*pi*j*m/n)
S[2*m+1] = sum_j=1^n/2 q[j]*sin(pi*j*(2*m+1)/n)
で計算すれば良い。そもそも、p, qを作っておかないと（つまりmが変わるごとに毎回a[j]+a[n-j]やa[j]-a[n-j]を計算していると）、[5]の変形による計算量の削減効果が発揮できません。

お礼

返答ありがとうございます。非常に丁寧に説明してくださいまして、とても感激しました。まだ完全に理解したわけではないので、質問を締め切ったのち、じっくり読ませていただきます。