§E20.16離散 Fourier 変換と FFT

最終更新

直交多項式による近似では、重み関数による積分で内積を定め、関数を直交系の係数に展開した。この積分を、NN個の等間隔点j/Nj/N（0≤j<N0\le j<N）に等しい重み1/N1/Nを置いた和に置き換えると、複素指数関数e2πikte^{2\pi\mathrm ikt}（0≤k<N0\le k<N）はこの和について正規直交系をなし、標本値xj:=f(j/N)x_j:=f(j/N)に対するe2πikte^{2\pi\mathrm ikt}の展開の係数はx^k/N\hat x_k/Nである。ここでx^k=∑j=0N−1xjωNjk\hat x_k=\sum_{j=0}^{N-1}x_j\omega_N^{jk}（ωN=e−2πi/N\omega_N=e^{-2\pi\mathrm i/N}）は長さNNの複素数列xxの離散 Fourier 変換であり、一の冪根の直交性によって、逆変換を明示的な式としてもつ。

離散 Fourier 変換は循環畳み込みを成分ごとの積に変える。多項式の積の係数を与える線形畳み込みも、長さmmとnnの二つの列を長さm+n−1m+n-1以上まで零で詰めれば、循環畳み込みとして計算することができる。定数ωNr\omega_N^rを計算済みとすると、離散 Fourier 変換を定義どおりに計算する方法は複素数の乗算をN2N^2回実行し、N=2sN=2^sのとき radix-2 高速 Fourier 変換（FFT）は同じ値を(N/2)log⁡2N(N/2)\log_2N回の乗算で計算する。

本記事では、離散 Fourier 変換と FFT について、基本的な性質と代表的な例を解説する。

1 一の冪根と離散 Fourier 変換

補題 1.1.N∈N≥1N\in\NNとし、ωN:=e−2πi/N\omega_N:=e^{-2\pi\mathrm i/N}と置く。任意のk∈Zk\in\Zに対して

∑j=0N−1ωNjk={N(N∣k),0(N∤k)\sum_{j=0}^{N-1}\omega_N^{jk}=\begin{cases}N&(N\mid k),\\0&(N\nmid k)\end{cases}

が成り立つ。

証明.k∈Zk\in\Zをとり、q:=ωNk=e−2πik/Nq:=\omega_N^k=e^{-2\pi\mathrm ik/N}と置く。N∣kN\mid kならばq=1q=1であり、和の各項は11であるから和はNNである。N∤kN\nmid kならばk/N∉Zk/N\notin\Zであり、eiθ=1e^{\mathrm i\theta}=1となる実数θ\thetaは2πZ2\pi\Zの元に限るからq≠1q\ne1である。qN=e−2πik=1q^N=e^{-2\pi\mathrm ik}=1であるから

∑j=0N−1qj=qN−1q−1=0\sum_{j=0}^{N-1}q^j=\frac{q^N-1}{q-1}=0

が成り立つ。▨

定義 1.2.N∈N≥1N\in\NNとし、ωN:=e−2πi/N\omega_N:=e^{-2\pi\mathrm i/N}と置く。CN\C^Nの元の成分を0,…,N−10,\dots,N-1で番号付ける。x=(x0,…,xN−1)∈CNx=(x_0,\dots,x_{N-1})\in\C^Nとk∈Zk\in\Zに対して

x^k:=∑j=0N−1xjωNjk\hat x_k:=\sum_{j=0}^{N-1}x_j\omega_N^{jk}

と置く。ωNN=1\omega_N^N=1であるから、任意のk∈Zk\in\Zに対してx^k+N=x^k\hat x_{k+N}=\hat x_kである。x^:=(x^0,…,x^N−1)∈CN\hat x:=(\hat x_0,\dots,\hat x_{N-1})\in\C^Nをxxの 離散 Fourier 変換 (discrete Fourier transform) といい、写像FN ⁣:CN→CN\mathcal F_N\colon\C^N\to\C^N、x↦x^x\mapsto\hat xを長さNNの離散 Fourier 変換という。

定理 1.3.N∈N≥1N\in\NNとし、CN\C^Nに標準内積⟨x,y⟩:=∑j=0N−1xjyj‾\langle x,y\rangle:=\sum_{j=0}^{N-1}x_j\overline{y_j}とノルム∥x∥2:=⟨x,x⟩1/2\lVert x\rVert_2:=\langle x,x\rangle^{1/2}を入れる。0≤k<N0\le k<Nに対してek∈CNe_k\in\C^Nを(ek)j:=N−1/2ωN−jk(e_k)_j:=N^{-1/2}\omega_N^{-jk}（0≤j<N0\le j<N）で定める。

  1. (e0,…,eN−1)(e_0,\dots,e_{N-1})はCN\C^Nの正規直交基底であり、任意のx∈CNx\in\C^Nと0≤k<N0\le k<Nに対して⟨x,ek⟩=N−1/2x^k\langle x,e_k\rangle=N^{-1/2}\hat x_kが成り立つ。
  2. 任意のx∈CNx\in\C^Nと0≤j<N0\le j<Nに対して xj=1N∑k=0N−1x^k ωN−jkx_j=\frac1N\sum_{k=0}^{N-1}\hat x_k\,\omega_N^{-jk} が成り立つ。FN\mathcal F_NはC\C上の線形同型である。
  3. 任意のx,y∈CNx,y\in\C^Nに対して⟨x^,y^⟩=N⟨x,y⟩\langle\hat x,\hat y\rangle=N\langle x,y\rangleと∥x^∥22=N∥x∥22\lVert\hat x\rVert_2^2=N\lVert x\rVert_2^2が成り立つ。

証明.0≤k,l<N0\le k,l<Nをとる。∣ωN∣=1|\omega_N|=1であるからωN−jl‾=ωNjl\overline{\omega_N^{-jl}}=\omega_N^{jl}であり、

⟨ek,el⟩=1N∑j=0N−1ωN−jkωNjl=1N∑j=0N−1ωNj(l−k)\langle e_k,e_l\rangle=\frac1N\sum_{j=0}^{N-1}\omega_N^{-jk}\omega_N^{jl}=\frac1N\sum_{j=0}^{N-1}\omega_N^{j(l-k)}

である。∣l−k∣<N|l-k|<NであるからN∣l−kN\mid l-kとk=lk=lは同値であり、補題 1.1により⟨ek,el⟩\langle e_k,e_l\rangleはk=lk=lのとき11、k≠lk\ne lのとき00である。c0,…,cN−1∈Cc_0,\dots,c_{N-1}\in\Cが∑kckek=0\sum_kc_ke_k=0を満たすならば、両辺とele_lの内積をとってcl=0c_l=0を得るので、e0,…,eN−1e_0,\dots,e_{N-1}は一次独立であり、dim⁡CCN=N\dim_\C\C^N=NであるからCN\C^Nの基底である。x∈CNx\in\C^Nに対して

⟨x,ek⟩=∑j=0N−1xjN−1/2ωN−jk‾=N−1/2∑j=0N−1xjωNjk=N−1/2x^k\langle x,e_k\rangle=\sum_{j=0}^{N-1}x_jN^{-1/2}\overline{\omega_N^{-jk}}=N^{-1/2}\sum_{j=0}^{N-1}x_j\omega_N^{jk}=N^{-1/2}\hat x_k

である。これで(1)は示された。

(1)と§E3.32 命題 2.3によりx=∑k=0N−1⟨x,ek⟩ek=∑k=0N−1N−1/2x^kekx=\sum_{k=0}^{N-1}\langle x,e_k\rangle e_k=\sum_{k=0}^{N-1}N^{-1/2}\hat x_ke_kであり、第jj成分をとって逆変換の式を得る。定義からFN\mathcal F_NはC\C-線形である。逆変換の式によりx^=0\hat x=0ならばx=0x=0であるからFN\mathcal F_Nは単射であり、有限次元線形空間CN\C^Nからそれ自身への単射線形写像であるから全単射である。

x=∑k⟨x,ek⟩ekx=\sum_k\langle x,e_k\rangle e_kとy=∑l⟨y,el⟩ely=\sum_l\langle y,e_l\rangle e_lを内積へ代入し、内積の第一変数についての線形性、第二変数についての共役線形性、(ek)(e_k)の正規直交性を用いると

⟨x,y⟩=∑k,l=0N−1⟨x,ek⟩⟨y,el⟩‾⟨ek,el⟩=∑k=0N−1⟨x,ek⟩⟨y,ek⟩‾=1N∑k=0N−1x^ky^k‾=1N⟨x^,y^⟩\langle x,y\rangle=\sum_{k,l=0}^{N-1}\langle x,e_k\rangle\overline{\langle y,e_l\rangle}\langle e_k,e_l\rangle=\sum_{k=0}^{N-1}\langle x,e_k\rangle\overline{\langle y,e_k\rangle}=\frac1N\sum_{k=0}^{N-1}\hat x_k\overline{\hat y_k}=\frac1N\langle\hat x,\hat y\rangle

である。y=xy=xとしてノルムの等式を得る。▨

注意 1.4.N∈N≥1N\in\NN、f ⁣:R→Cf\colon\R\to\Cとし、xj:=f(j/N)x_j:=f(j/N)（0≤j<N0\le j<N）、φk(t):=e2πikt\varphi_k(t):=e^{2\pi\mathrm ikt}（0≤k<N0\le k<N）と置く。φk(j/N)=ωN−jk\varphi_k(j/N)=\omega_N^{-jk}であるから、定理 1.3 (2)により三角多項式p:=N−1∑k=0N−1x^kφkp:=N^{-1}\sum_{k=0}^{N-1}\hat x_k\varphi_kはp(j/N)=f(j/N)p(j/N)=f(j/N)（0≤j<N0\le j<N）を満たす。関数g,h ⁣:R→Cg,h\colon\R\to\Cに対して⟨g,h⟩N:=N−1∑j=0N−1g(j/N)h(j/N)‾\langle g,h\rangle_N:=N^{-1}\sum_{j=0}^{N-1}g(j/N)\overline{h(j/N)}と置く。(φk(j/N))0≤j<N=N1/2ek(\varphi_k(j/N))_{0\le j<N}=N^{1/2}e_kであるから、定理 1.3 (1)により⟨φk,φl⟩N=⟨ek,el⟩\langle\varphi_k,\varphi_l\rangle_N=\langle e_k,e_l\rangleはk=lk=lのとき11、k≠lk\ne lのとき00であり、⟨f,φk⟩N=x^k/N\langle f,\varphi_k\rangle_N=\hat x_k/Nはppにおけるφk\varphi_kの係数に等しい。g=∑kckφkg=\sum_kc_k\varphi_kが⟨g,g⟩N=0\langle g,g\rangle_N=0を満たすならば∑kckN1/2ek=0\sum_kc_kN^{1/2}e_k=0でありc=0c=0であるから、⟨ , ⟩N\langle\ ,\ \rangle_Nはφ0,…,φN−1\varphi_0,\dots,\varphi_{N-1}の張る空間の上で内積である。重み関数による積分で定める内積に対する多項式の直交展開と比べると、⟨ , ⟩N\langle\ ,\ \rangle_Nは積分をNN個の点j/Nj/Nに等しい重み1/N1/Nを置いた和に置き換えたものである。

2 循環畳み込み

定義 2.1.N∈N≥1N\in\NNとし、整数jjをNNで割った余りをj mod N∈{0,…,N−1}j\bmod N\in\{0,\dots,N-1\}と書く。x,y∈CNx,y\in\C^Nに対して、x⊛y∈CNx\circledast y\in\C^Nを

(x⊛y)j:=∑l=0N−1xl y(j−l) mod N(0≤j<N)(x\circledast y)_j:=\sum_{l=0}^{N-1}x_l\,y_{(j-l)\bmod N}\qquad(0\le j<N)

で定め、xxとyyの 循環畳み込み (circular convolution) という。

補題 2.2.N∈N≥1N\in\NN、x,y∈CNx,y\in\C^Nとし、0≤j<N0\le j<Nに対してSj:={(l,i)∈{0,…,N−1}2∣l+i≡j(modN)}S_j:=\{(l,i)\in\{0,\dots,N-1\}^2\mid l+i\equiv j\pmod N\}と置く。このとき

(x⊛y)j=∑(l,i)∈Sjxlyi(0≤j<N)(x\circledast y)_j=\sum_{(l,i)\in S_j}x_ly_i\qquad(0\le j<N)

が成り立つ。

証明. 演習とする（問題 6.1）。▨

定理 2.3.N∈N≥1N\in\NN、x,y∈CNx,y\in\C^Nとし、x⊙y:=(x0y0,…,xN−1yN−1)x\odot y:=(x_0y_0,\dots,x_{N-1}y_{N-1})と置く。このときx⊛y^=x^⊙y^\widehat{x\circledast y}=\hat x\odot\hat yであり、x⊛y=FN−1(x^⊙y^)x\circledast y=\mathcal F_N^{-1}(\hat x\odot\hat y)が成り立つ。

証明.0≤k<N0\le k<Nをとり、SjS_jを補題 2.2の集合とする。補題 2.2により

x⊛y^k=∑j=0N−1∑(l,i)∈Sjxlyi ωNjk\widehat{x\circledast y}_k=\sum_{j=0}^{N-1}\sum_{(l,i)\in S_j}x_ly_i\,\omega_N^{jk}

である。(l,i)∈Sj(l,i)\in S_jならばある整数rrについてl+i=j+rNl+i=j+rNであり、ωNN=1\omega_N^N=1であるからωNjk=ωN(l+i)k\omega_N^{jk}=\omega_N^{(l+i)k}である。各(l,i)∈{0,…,N−1}2(l,i)\in\{0,\dots,N-1\}^2はj=(l+i) mod Nj=(l+i)\bmod Nに対するSjS_jだけに属するので、S0,…,SN−1S_0,\dots,S_{N-1}は{0,…,N−1}2\{0,\dots,N-1\}^2の分割である。したがって

x⊛y^k=∑l=0N−1∑i=0N−1xlyi ωNlkωNik=x^ky^k\widehat{x\circledast y}_k=\sum_{l=0}^{N-1}\sum_{i=0}^{N-1}x_ly_i\,\omega_N^{lk}\omega_N^{ik}=\hat x_k\hat y_k

である。定理 1.3 (2)によりFN\mathcal F_Nは全単射であるから、x⊛y=FN−1(x^⊙y^)x\circledast y=\mathcal F_N^{-1}(\hat x\odot\hat y)である。▨

3 高速 Fourier 変換

補題 3.1.M∈N≥1M\in\NN、N:=2MN:=2M、x∈CNx\in\C^Nとし、a,b∈CMa,b\in\C^Mをam:=x2ma_m:=x_{2m}、bm:=x2m+1b_m:=x_{2m+1}（0≤m<M0\le m<M）で定める。a^,b^\hat a,\hat bを長さMMの離散 Fourier 変換とする。このとき0≤k<M0\le k<Mに対して

x^k=a^k+ωNkb^k,x^k+M=a^k−ωNkb^k\hat x_k=\hat a_k+\omega_N^k\hat b_k,\qquad\hat x_{k+M}=\hat a_k-\omega_N^k\hat b_k

が成り立つ。

証明.κ∈Z\kappa\in\Zをとる。和をj=2mj=2mとj=2m+1j=2m+1（0≤m<M0\le m<M）に分けると

x^κ=∑m=0M−1x2mωN2mκ+ωNκ∑m=0M−1x2m+1ωN2mκ\hat x_\kappa=\sum_{m=0}^{M-1}x_{2m}\omega_N^{2m\kappa}+\omega_N^\kappa\sum_{m=0}^{M-1}x_{2m+1}\omega_N^{2m\kappa}

である。ωN2=e−2πi/M=ωM\omega_N^2=e^{-2\pi\mathrm i/M}=\omega_Mであるから、右辺はa^κ+ωNκb^κ\hat a_\kappa+\omega_N^\kappa\hat b_\kappaに等しい。κ=k\kappa=kとして第一式を得る。κ=k+M\kappa=k+Mとすると、定義 1.2によりa^k+M=a^k\hat a_{k+M}=\hat a_k、b^k+M=b^k\hat b_{k+M}=\hat b_kであり、ωNM=e−πi=−1\omega_N^M=e^{-\pi\mathrm i}=-1であるからωNk+M=−ωNk\omega_N^{k+M}=-\omega_N^kであり、第二式を得る。▨

定義 3.2.s∈N≥0s\in\Nとし、N:=2sN:=2^sと置く。複素数ωLk\omega_L^k（L=2rL=2^r、1≤r≤s1\le r\le s、0≤k<L/20\le k<L/2）の値は計算済みの定数として与えられるものとする。写像ΦN ⁣:CN→CN\Phi_N\colon\C^N\to\C^Nと、ΦN(x)\Phi_N(x)の計算で実行する複素数の演算を、ssに関して帰納的に次のように定める。

  1. s=0s=0のとき、Φ1(x):=x\Phi_1(x):=xとし、演算を実行しない。
  2. s≥1s\ge1のとき、M:=N/2M:=N/2と置き、a,b∈CMa,b\in\C^Mをam:=x2ma_m:=x_{2m}、bm:=x2m+1b_m:=x_{2m+1}で定め、α:=ΦM(a)\alpha:=\Phi_M(a)、β:=ΦM(b)\beta:=\Phi_M(b)を計算する。0≤k<M0\le k<Mの各kkについて、乗算tk:=ωNkβkt_k:=\omega_N^k\beta_kを一回、加算αk+tk\alpha_k+t_kと減算αk−tk\alpha_k-t_kを一回ずつ実行し、ΦN(x)k:=αk+tk\Phi_N(x)_k:=\alpha_k+t_k、ΦN(x)k+M:=αk−tk\Phi_N(x)_{k+M}:=\alpha_k-t_kと置く。ΦN(x)\Phi_N(x)の計算の演算は、α\alphaとβ\betaの計算の演算と、これらの3M3M回の演算である。

ΦN\Phi_Nを長さNNの radix-2 高速 Fourier 変換 (radix-2 fast Fourier transform)（radix-2 FFT）という。

定理 3.3.s∈N≥0s\in\Nとし、N:=2sN:=2^sと置く。

  1. 任意のx∈CNx\in\C^Nに対してΦN(x)=x^\Phi_N(x)=\hat xである。
  2. ΦN(x)\Phi_N(x)の計算は、複素数の乗算をちょうど(N/2)log⁡2N(N/2)\log_2N回、複素数の加算・減算をちょうどNlog⁡2NN\log_2N回実行する。複素数の乗算を実数の乗算44回と加減算22回で、複素数の加減算を実数の加減算22回で実行するとき、実数の演算の総数は5Nlog⁡2N5N\log_2Nである。
  3. 定数ωNr\omega_N^r（0≤r<N0\le r<N）を与え、0≤k<N0\le k<Nの各kkについてx^k\hat x_kをNN個の積xjωN(jk) mod Nx_j\omega_N^{(jk)\bmod N}のN−1N-1回の加算による和として計算すると、複素数の乗算はN2N^2回、加算はN(N−1)N(N-1)回であり、前項と同じ換算で実数の演算の総数は8N2−2N8N^2-2Nである。

証明.(1)を示す。s=0s=0ならば、任意のx∈C1x\in\C^1に対してx^0=x0ω10=x0=Φ1(x)0\hat x_0=x_0\omega_1^0=x_0=\Phi_1(x)_0である。s≥1s\ge1とし、M:=N/2M:=N/2についてΦM=FM\Phi_M=\mathcal F_Mが成り立つとする。x∈CNx\in\C^Nをとり、a,ba,bを定義 3.2のとおりにとるとα=a^\alpha=\hat a、β=b^\beta=\hat bであり、0≤k<M0\le k<Mに対してΦN(x)k=a^k+ωNkb^k\Phi_N(x)_k=\hat a_k+\omega_N^k\hat b_k、ΦN(x)k+M=a^k−ωNkb^k\Phi_N(x)_{k+M}=\hat a_k-\omega_N^k\hat b_kである。補題 3.1により右辺はそれぞれx^k\hat x_k、x^k+M\hat x_{k+M}に等しい。{k,k+M∣0≤k<M}={0,…,N−1}\{k,k+M\mid 0\le k<M\}=\{0,\dots,N-1\}であるからΦN(x)=x^\Phi_N(x)=\hat xであり、ssに関する帰納法により任意のs∈N≥0s\in\NについてΦ2s=F2s\Phi_{2^s}=\mathcal F_{2^s}である。

(2)を示す。Φ2s\Phi_{2^s}の乗算の回数をμs\mu_s、加減算の回数をσs\sigma_sと書く。定義 3.2によりμ0=σ0=0\mu_0=\sigma_0=0であり、s≥1s\ge1に対して

μs=2μs−1+2s−1,σs=2σs−1+2s\mu_s=2\mu_{s-1}+2^{s-1},\qquad\sigma_s=2\sigma_{s-1}+2^s

である。μs−1=(s−1)2s−2\mu_{s-1}=(s-1)2^{s-2}、σs−1=(s−1)2s−1\sigma_{s-1}=(s-1)2^{s-1}ならばμs=(s−1)2s−1+2s−1=s2s−1\mu_s=(s-1)2^{s-1}+2^{s-1}=s2^{s-1}、σs=(s−1)2s+2s=s2s\sigma_s=(s-1)2^s+2^s=s2^sであり、ssに関する帰納法によりμs=(N/2)s\mu_s=(N/2)s、σs=Ns\sigma_s=Nsである。実数の演算の総数は6μs+2σs=3Ns+2Ns=5Ns6\mu_s+2\sigma_s=3Ns+2Ns=5Nsである。

(3)を示す。各kkについて乗算はNN回、加算はN−1N-1回であり、kkはNN通りである。実数の演算の総数は6N2+2N(N−1)=8N2−2N6N^2+2N(N-1)=8N^2-2Nである。▨

注意 3.4.s∈N≥1s\in\NN、N:=2sN:=2^sとする。ΦN\Phi_Nの再帰で長さL=2rL=2^r（1≤r≤s1\le r\le s）の段が用いる定数ωLk\omega_L^k（0≤k<L/20\le k<L/2）は、ωL=e−2πi/L=ωNN/L\omega_L=e^{-2\pi\mathrm i/L}=\omega_N^{N/L}によりωNkN/L\omega_N^{kN/L}に等しく、0≤kN/L<N/20\le kN/L<N/2である。したがってすべての段の定数はωN0,…,ωNN/2−1\omega_N^0,\dots,\omega_N^{N/2-1}のN/2N/2個から取ることができる。定理 3.3 (2)の回数はこれらの定数の計算を含まない。

例 3.5.N=4N=4、x=(1,2,3,0)x=(1,2,3,0)とする。ω4=e−πi/2=−i\omega_4=e^{-\pi\mathrm i/2}=-\mathrm i、ω2=−1\omega_2=-1である。a=(x0,x2)=(1,3)a=(x_0,x_2)=(1,3)、b=(x1,x3)=(2,0)b=(x_1,x_3)=(2,0)であり、長さ22の段でα=Φ2(a)=(1+3,1−3)=(4,−2)\alpha=\Phi_2(a)=(1+3,1-3)=(4,-2)、β=Φ2(b)=(2,2)\beta=\Phi_2(b)=(2,2)である。長さ44の段ではt0=ω40β0=2t_0=\omega_4^0\beta_0=2、t1=ω4β1=−2it_1=\omega_4\beta_1=-2\mathrm iであり、

Φ4(x)=(α0+t0, α1+t1, α0−t0, α1−t1)=(6, −2−2i, 2, −2+2i)\Phi_4(x)=(\alpha_0+t_0,\ \alpha_1+t_1,\ \alpha_0-t_0,\ \alpha_1-t_1)=(6,\ -2-2\mathrm i,\ 2,\ -2+2\mathrm i)

である。定義による計算x^1=1+2(−i)+3(−i)2=−2−2i\hat x_1=1+2(-\mathrm i)+3(-\mathrm i)^2=-2-2\mathrm i、x^2=1−2+3=2\hat x_2=1-2+3=2、x^3=1+2i−3=−2+2i\hat x_3=1+2\mathrm i-3=-2+2\mathrm i、x^0=6\hat x_0=6と一致する。乗算は長さ22の段で22回、長さ44の段で22回の計4=(4/2)log⁡244=(4/2)\log_24回であり、加減算は4+4=8=4log⁡244+4=8=4\log_24回である。∥x∥22=1+4+9=14\lVert x\rVert_2^2=1+4+9=14、∥x^∥22=36+8+4+8=56=4⋅14\lVert\hat x\rVert_2^2=36+8+4+8=56=4\cdot14であり、定理 1.3 (3)の等式が成り立っている。

注意 3.6.N∈N≥1N\in\NN、x∈CN∖{0}x\in\C^N\setminus\{0\}、h∈CNh\in\C^Nとする。定義 1.2によりFN\mathcal F_Nは線形であるからFN(x+h)−FNx=FNh\mathcal F_N(x+h)-\mathcal F_Nx=\mathcal F_Nhであり、定理 1.3 (3)により∥FNh∥2=N1/2∥h∥2\lVert\mathcal F_Nh\rVert_2=N^{1/2}\lVert h\rVert_2、∥FNx∥2=N1/2∥x∥2\lVert\mathcal F_Nx\rVert_2=N^{1/2}\lVert x\rVert_2である。x≠0x\ne0であるから∥FNx∥2>0\lVert\mathcal F_Nx\rVert_2>0であり、

∥FN(x+h)−FNx∥2∥FNx∥2=N1/2∥h∥2N1/2∥x∥2=∥h∥2∥x∥2\frac{\lVert\mathcal F_N(x+h)-\mathcal F_Nx\rVert_2}{\lVert\mathcal F_Nx\rVert_2}=\frac{N^{1/2}\lVert h\rVert_2}{N^{1/2}\lVert x\rVert_2}=\frac{\lVert h\rVert_2}{\lVert x\rVert_2}

が成り立つ。

例 3.7. CPython の float（基数22、仮数5353桁、最近接偶数丸め、u=2−53≈1.11×10−16u=2^{-53}\approx1.11\times10^{-16}）で計算する。s∈{6,8,10,12}s\in\{6,8,10,12\}、N=2sN=2^sとし、整数の入力xj:=((37j2+11j) mod 101)−50x_j:=((37j^2+11j)\bmod101)-50（0≤j<N0\le j<N）をとる。

定数ωNr\omega_N^r（0≤r<N0\le r<N）は、実部を math.cos(2*math.pi*r/N)、虚部を -math.sin(2*math.pi*r/N) で計算した値ω~Nr\tilde\omega_N^rに置き換え、複素数の積(p+qi)(p′+q′i)(p+q\mathrm i)(p'+q'\mathrm i)は(pp′−qq′)+(pq′+qp′)i(pp'-qq')+(pq'+qp')\mathrm iの各演算を丸めて計算する。直接計算は、各kkについて積xjω~N(jk) mod Nx_j\tilde\omega_N^{(jk)\bmod N}をj=0,1,…,N−1j=0,1,\dots,N-1の順に逐次和で加え、FFT はΦN\Phi_NのωLk\omega_L^kを注意 3.4のω~NkN/L\tilde\omega_N^{kN/L}に置き換えて計算する。

計算値y~\tilde yと、5050桁の十進演算で計算したx^\hat xから相対誤差∥y~−x^∥2/∥x^∥2\lVert\tilde y-\hat x\rVert_2/\lVert\hat x\rVert_2を求めると次のとおりである。

NN 直接計算 FFT
262^6 3.1×10−163.1\times10^{-16} 2.8×10−162.8\times10^{-16}
282^8 6.1×10−166.1\times10^{-16} 3.2×10−163.2\times10^{-16}
2102^{10} 1.0×10−151.0\times10^{-15} 3.9×10−163.9\times10^{-16}
2122^{12} 2.1×10−152.1\times10^{-15} 4.7×10−164.7\times10^{-16}

表の値は一つの入力に対する観察であり、誤差の上界ではない。

4 零詰めと線形畳み込み

定義 4.1.m,n∈N≥1m,n\in\NN、x∈Cmx\in\C^m、y∈Cny\in\C^nとする。0≤p≤m+n−20\le p\le m+n-2に対してTp:={(l,i)∣0≤l<m, 0≤i<n, l+i=p}T_p:=\{(l,i)\mid 0\le l<m,\ 0\le i<n,\ l+i=p\}と置き、x∗y∈Cm+n−1x*y\in\C^{m+n-1}を(x∗y)p:=∑(l,i)∈Tpxlyi(x*y)_p:=\sum_{(l,i)\in T_p}x_ly_iで定めて、xxとyyの 線形畳み込み (linear convolution) という。整数N≥mN\ge mに対して、x[N]:=(x0,…,xm−1,0,…,0)∈CNx^{[N]}:=(x_0,\dots,x_{m-1},0,\dots,0)\in\C^Nをxxの長さNNへの 零詰め (zero padding) という。

命題 4.2.m,n,N∈N≥1m,n,N\in\NNとし、N≥max⁡(m,n)N\ge\max(m,n)、x∈Cmx\in\C^m、y∈Cny\in\C^n、z:=x∗yz:=x*yとする。

  1. 0≤j<N0\le j<Nに対して (x[N]⊛y[N])j=∑0≤p≤m+n−2p≡j (mod N)zp(x^{[N]}\circledast y^{[N]})_j=\sum_{\substack{0\le p\le m+n-2\\ p\equiv j\ (\mathrm{mod}\ N)}}z_p が成り立つ。
  2. N≥m+n−1N\ge m+n-1ならば、0≤j≤m+n−20\le j\le m+n-2に対して(x[N]⊛y[N])j=zj(x^{[N]}\circledast y^{[N]})_j=z_jであり、m+n−1≤j<Nm+n-1\le j<Nに対して(x[N]⊛y[N])j=0(x^{[N]}\circledast y^{[N]})_j=0である。

証明.0≤j<N0\le j<Nをとり、x′:=x[N]x':=x^{[N]}、y′:=y[N]y':=y^{[N]}とする。補題 2.2により(x′⊛y′)j=∑(l,i)∈Sjxl′yi′(x'\circledast y')_j=\sum_{(l,i)\in S_j}x'_ly'_iである。l≥ml\ge mまたはi≥ni\ge nである項は00であるから、和は0≤l<m0\le l<m、0≤i<n0\le i<n、l+i≡j(modN)l+i\equiv j\pmod Nを満たす(l,i)(l,i)の上の和∑xlyi\sum x_ly_iに等しい。この集合は、0≤p≤m+n−20\le p\le m+n-2かつp≡j(modN)p\equiv j\pmod Nを満たすppに対するTpT_pの交わらない和集合であり、(1)を得る。N≥m+n−1N\ge m+n-1ならば、(1)の和に現れるppは0≤p<N0\le p<Nを満たし、0≤j<N0\le j<Nとp≡j(modN)p\equiv j\pmod Nからp=jp=jである。j≤m+n−2j\le m+n-2ならば和はzjz_jであり、j≥m+n−1j\ge m+n-1ならば和は空であるから00である。▨

系 4.3.m,n∈N≥1m,n\in\NN、s∈N≥0s\in\N、N:=2s≥m+n−1N:=2^s\ge m+n-1、x∈Cmx\in\C^m、y∈Cny\in\C^nとする。w:=ΦN(x[N])⊙ΦN(y[N])w:=\Phi_N(x^{[N]})\odot\Phi_N(y^{[N]})と置き、w‾\overline wを成分ごとの複素共役とする。このとき0≤p≤m+n−20\le p\le m+n-2に対して

(x∗y)p=1N ΦN(w‾)p‾(x*y)_p=\frac1N\,\overline{\Phi_N(\overline w)_p}

が成り立つ。ΦN\Phi_Nの三回の計算、⊙\odotのNN回の乗算、2N2N回の複素共役、実数1/N1/NによるNN回の乗算によって右辺を0≤p<N0\le p<Nについて計算すると、複素数の乗算は(3/2)Nlog⁡2N+N(3/2)N\log_2N+N回、加減算は3Nlog⁡2N3N\log_2N回である。ssを2s≥m+n−12^s\ge m+n-1を満たす最小の整数にとるとN<2(m+n−1)N<2(m+n-1)である。定義によるx∗yx*yの計算は、乗算をmnmn回、加算をmn−(m+n−1)mn-(m+n-1)回実行する。

証明.x′:=x[N]x':=x^{[N]}、y′:=y[N]y':=y^{[N]}、v:=x′⊛y′v:=x'\circledast y'と置く。定理 3.3 (1)と定理 2.3によりw=x′^⊙y′^=v^w=\hat{x'}\odot\hat{y'}=\hat vである。定理 1.3 (2)により、0≤p<N0\le p<Nに対して

vp=1N∑k=0N−1wkωN−pk=1N ∑k=0N−1wk‾ ωNpk‾=1N w‾^p‾v_p=\frac1N\sum_{k=0}^{N-1}w_k\omega_N^{-pk}=\frac1N\,\overline{\sum_{k=0}^{N-1}\overline{w_k}\,\omega_N^{pk}}=\frac1N\,\overline{\widehat{\overline w}_p}

であり、定理 3.3 (1)によりw‾^=ΦN(w‾)\widehat{\overline w}=\Phi_N(\overline w)である。N≥m+n−1N\ge m+n-1であるから、命題 4.2 (2)により0≤p≤m+n−20\le p\le m+n-2に対してvp=(x∗y)pv_p=(x*y)_pである。演算の回数は定理 3.3 (2)を三回分加え、⊙\odotのNN回を加えたものである。最小のssについて、s=0s=0ならばN=1<2≤2(m+n−1)N=1<2\le2(m+n-1)であり、s≥1s\ge1ならば2s−1<m+n−12^{s-1}<m+n-1であるからN<2(m+n−1)N<2(m+n-1)である。定義による計算では、各ppについて∣Tp∣|T_p|回の乗算と∣Tp∣−1|T_p|-1回の加算を実行する。0≤p≤m+n−20\le p\le m+n-2に対して(max⁡(0,p−n+1), p−max⁡(0,p−n+1))∈Tp(\max(0,p-n+1),\,p-\max(0,p-n+1))\in T_pであるからTp≠∅T_p\ne\emptysetであり、T0,…,Tm+n−2T_0,\dots,T_{m+n-2}は{0,…,m−1}×{0,…,n−1}\{0,\dots,m-1\}\times\{0,\dots,n-1\}の分割であるから∑p∣Tp∣=mn\sum_p|T_p|=mnである。加算の回数はmn−(m+n−1)mn-(m+n-1)である。▨

例 4.4.m=3m=3、n=2n=2、x=(1,2,3)x=(1,2,3)、y=(1,1)y=(1,1)とする。z=x∗y=(1,3,5,3)z=x*y=(1,3,5,3)であり、これは多項式の積(1+2t+3t2)(1+t)=1+3t+5t2+3t3(1+2t+3t^2)(1+t)=1+3t+5t^2+3t^3の係数である。N=4=m+n−1N=4=m+n-1とする。例 3.5によりx[4]^=(6,−2−2i,2,−2+2i)\widehat{x^{[4]}}=(6,-2-2\mathrm i,2,-2+2\mathrm i)であり、y[4]^=(1+1, 1−i, 1−1, 1+i)=(2,1−i,0,1+i)\widehat{y^{[4]}}=(1+1,\ 1-\mathrm i,\ 1-1,\ 1+\mathrm i)=(2,1-\mathrm i,0,1+\mathrm i)である。成分ごとの積はw=(12,−4,0,−4)w=(12,-4,0,-4)であり、ω4−1=i\omega_4^{-1}=\mathrm iであるから定理 2.3と定理 1.3 (2)により

(x[4]⊛y[4])j=14(12−4 i j−4 i 3j)=3−i j−i 3j(0≤j<4)(x^{[4]}\circledast y^{[4]})_j=\frac14\bigl(12-4\,\mathrm i^{\,j}-4\,\mathrm i^{\,3j}\bigr)=3-\mathrm i^{\,j}-\mathrm i^{\,3j}\qquad(0\le j<4)

である。j=0,1,2,3j=0,1,2,3で値は1,3,5,31,3,5,3であり、zzに一致する。N=3N=3とするとN≥max⁡(3,2)N\ge\max(3,2)であるがN<m+n−1N<m+n-1であり、定義からx[3]⊛y[3]=(1⋅1+3⋅1, 1⋅1+2⋅1, 2⋅1+3⋅1)=(4,3,5)x^{[3]}\circledast y^{[3]}=(1\cdot1+3\cdot1,\ 1\cdot1+2\cdot1,\ 2\cdot1+3\cdot1)=(4,3,5)である。これは命題 4.2 (1)の(z0+z3,z1,z2)(z_0+z_3,z_1,z_2)であり、t3t^3の係数z3z_3が添字00へ折り返されている。

5 区別することができない周波数

命題 5.1.N∈N≥1N\in\NN、ν,ν′∈R\nu,\nu'\in\Rとする。条件PPに対して、[P][P]をPPが成り立つとき11、成り立たないとき00と置く。

  1. 任意のj∈Zj\in\Zに対してe2πiνj/N=e2πiν′j/Ne^{2\pi\mathrm i\nu j/N}=e^{2\pi\mathrm i\nu'j/N}が成り立つことと、ν′−ν∈NZ\nu'-\nu\in N\Zであることは同値である。
  2. 任意のj∈Zj\in\Zに対してcos⁡(2πνj/N)=cos⁡(2πν′j/N)\cos(2\pi\nu j/N)=\cos(2\pi\nu'j/N)が成り立つことと、ν′−ν∈NZ\nu'-\nu\in N\Zまたはν′+ν∈NZ\nu'+\nu\in N\Zであることは同値である。
  3. ν∈Z\nu\in\Zとし、x,c∈CNx,c\in\C^Nをxj:=e2πiνj/Nx_j:=e^{2\pi\mathrm i\nu j/N}、cj:=cos⁡(2πνj/N)c_j:=\cos(2\pi\nu j/N)（0≤j<N0\le j<N）で定める。このとき0≤k<N0\le k<Nに対して x^k=N [ k≡ν (mod N) ],c^k=N2([ k≡ν (mod N) ]+[ k≡−ν (mod N) ])\hat x_k=N\,[\,k\equiv\nu\ (\mathrm{mod}\ N)\,],\qquad\hat c_k=\frac N2\bigl([\,k\equiv\nu\ (\mathrm{mod}\ N)\,]+[\,k\equiv-\nu\ (\mathrm{mod}\ N)\,]\bigr) が成り立つ。

証明.(1)を示す。ν′−ν=rN\nu'-\nu=rN（r∈Zr\in\Z）ならば、任意のj∈Zj\in\Zに対してe2πiν′j/N=e2πiνj/Ne2πirj=e2πiνj/Ne^{2\pi\mathrm i\nu'j/N}=e^{2\pi\mathrm i\nu j/N}e^{2\pi\mathrm irj}=e^{2\pi\mathrm i\nu j/N}である。逆にj=1j=1で等式が成り立つならばe2πi(ν′−ν)/N=1e^{2\pi\mathrm i(\nu'-\nu)/N}=1であり、eiθ=1e^{\mathrm i\theta}=1となる実数θ\thetaは2πZ2\pi\Zの元に限るから(ν′−ν)/N∈Z(\nu'-\nu)/N\in\Zである。

(2)を示す。実数θ,θ′\theta,\theta'についてcos⁡θ−cos⁡θ′=−2sin⁡θ+θ′2sin⁡θ−θ′2\cos\theta-\cos\theta'=-2\sin\frac{\theta+\theta'}2\sin\frac{\theta-\theta'}2であり、sin⁡\sinの零点はπZ\pi\Zであるから、cos⁡θ=cos⁡θ′\cos\theta=\cos\theta'と「θ′−θ∈2πZ\theta'-\theta\in2\pi\Zまたはθ′+θ∈2πZ\theta'+\theta\in2\pi\Z」は同値である。ν′=±ν+rN\nu'=\pm\nu+rN（r∈Zr\in\Z）ならば、任意のj∈Zj\in\Zに対して2πν′j/N=±2πνj/N+2πrj2\pi\nu'j/N=\pm2\pi\nu j/N+2\pi rjであり、余弦の値は等しい。逆にj=1j=1で等式が成り立つならば、θ:=2πν/N\theta:=2\pi\nu/N、θ′:=2πν′/N\theta':=2\pi\nu'/Nに上の同値を適用して、ν′−ν∈NZ\nu'-\nu\in N\Zまたはν′+ν∈NZ\nu'+\nu\in N\Zを得る。

(3)を示す。xj=ωN−νjx_j=\omega_N^{-\nu j}であるからx^k=∑j=0N−1ωNj(k−ν)\hat x_k=\sum_{j=0}^{N-1}\omega_N^{j(k-\nu)}であり、補題 1.1により第一式を得る。xj′:=e−2πiνj/Nx'_j:=e^{-2\pi\mathrm i\nu j/N}と置くとc=(x+x′)/2c=(x+x')/2であり、第一式を−ν-\nuに適用したx′^k=N[ k≡−ν (mod N) ]\hat{x'}_k=N[\,k\equiv-\nu\ (\mathrm{mod}\ N)\,]とFN\mathcal F_Nの線形性から第二式を得る。▨

例 5.2.N=8N=8、ν=3\nu=3とする。ν′=5\nu'=5ではν′+ν=8∈8Z\nu'+\nu=8\in8\Z、ν′−ν=2∉8Z\nu'-\nu=2\notin8\Zであるから、命題 5.1 (2)により余弦の標本は一致し、命題 5.1 (1)により複素指数の標本は一致しない。実際j=1j=1でe2πi⋅3/8=e3πi/4≠e5πi/4=e2πi⋅5/8e^{2\pi\mathrm i\cdot3/8}=e^{3\pi\mathrm i/4}\ne e^{5\pi\mathrm i/4}=e^{2\pi\mathrm i\cdot5/8}である。ν′=11\nu'=11ではν′−ν=8\nu'-\nu=8であり、複素指数の標本も余弦の標本も一致する。命題 5.1 (3)により、ν=3\nu=3の複素指数の標本の離散 Fourier 変換はk=3k=3だけで値88をとり、ν=5\nu=5のものはk=5k=5だけで値88をとる。余弦の標本の離散 Fourier 変換は、ν=3,5,11,13\nu=3,5,11,13のいずれについてもk=3k=3とk=5k=5で値44、他のkkで00である。

6 演習

問題 6.1.補題 2.2の証明を完成させよ。

解答.

0≤j<N0\le j<Nをとる。l,i∈{0,…,N−1}l,i\in\{0,\dots,N-1\}について、l+i≡j(modN)l+i\equiv j\pmod Nはi≡j−l(modN)i\equiv j-l\pmod Nと同値であり、{0,…,N−1}\{0,\dots,N-1\}の中でj−lj-lとNNを法として合同な元は(j−l) mod N(j-l)\bmod Nだけであるから、(l,i)∈Sj(l,i)\in S_jとi=(j−l) mod Ni=(j-l)\bmod Nは同値である。したがってl↦(l,(j−l) mod N)l\mapsto(l,(j-l)\bmod N)は{0,…,N−1}\{0,\dots,N-1\}からSjS_jへの全単射であり、逆写像は(l,i)↦l(l,i)\mapsto lである。この全単射で和の添字を付け替えると

∑(l,i)∈Sjxlyi=∑l=0N−1xl y(j−l) mod N=(x⊛y)j\sum_{(l,i)\in S_j}x_ly_i=\sum_{l=0}^{N-1}x_l\,y_{(j-l)\bmod N}=(x\circledast y)_j

である。▨

問題 6.2.s∈N≥0s\in\N、N:=2sN:=2^sとする。定義 3.2のΦN\Phi_Nの計算で各段が実行する乗算tkt_kのうち、k≠0k\ne0であるものの回数は(N/2)log⁡2N−N+1(N/2)\log_2N-N+1であることを示せ。

解答.

Φ2s\Phi_{2^s}の計算でk≠0k\ne0である乗算の回数をνs\nu_sと書く。ν0=0\nu_0=0である。s≥1s\ge1のとき、長さ2s2^sの段はk=0,…,2s−1−1k=0,\dots,2^{s-1}-1の2s−12^{s-1}回の乗算を実行し、そのうちk≠0k\ne0であるものは2s−1−12^{s-1}-1回であるから、νs=2νs−1+2s−1−1\nu_s=2\nu_{s-1}+2^{s-1}-1である。νs−1=(s−1)2s−2−2s−1+1\nu_{s-1}=(s-1)2^{s-2}-2^{s-1}+1ならば

νs=(s−1)2s−1−2s+2+2s−1−1=s2s−1−2s+1\nu_s=(s-1)2^{s-1}-2^s+2+2^{s-1}-1=s2^{s-1}-2^s+1

であり、s=0s=0でs2s−1−2s+1=0=ν0s2^{s-1}-2^s+1=0=\nu_0であるから、ssに関する帰納法によりνs=(N/2)log⁡2N−N+1\nu_s=(N/2)\log_2N-N+1である。▨

前提記事