#この記事の概要
FFTの原理をなるべく数式で書いて、Juliaで実装します。

また、この記事を書くに当たって以下の本、記事をめちゃめちゃ参考にしました。

物理数学の直観的方法 (https://bookclub.kodansha.co.jp/product?item=0000194699 )
EMANの物理学（https://eman-physics.net/math/fourier01.html ）
安居院猛, 中嶋正之, 『FFT の使い方』, 産報出版株式会社 (1982)
樋口龍雄, 川又政征, 『ディジタル信号処理-MATLAB 対応-』, 昭晃堂 (2000)
FFT(高速フーリエ変換):実装のための数式と実装例 (https://qiita.com/k-kotera/items/4fe20c48be56d2eb4ba4)

#DFTの動作
N個のデータが以下のようにベクトルの形で与えられるとする。

\boldsymbol{f} = 
    \begin{pmatrix}
      f_0 \\
      \vdots \\
      f_{N-1} \\
    \end{pmatrix}

この時、$\boldsymbol{f}$は直交基底ベクトルの組{$\boldsymbol{e}_0,\boldsymbol{e}_1,\boldsymbol{e}_2,...$}と、その係数の組{$f_0,f_1,f_2,...$}で以下のように表せる。

\boldsymbol{f} = f_0\boldsymbol{e}_0+f_1\boldsymbol{e}_1+...+f_{N-1}\boldsymbol{e}_{N-1}

ただし、{$\boldsymbol{e}_0,\boldsymbol{e}_1,\boldsymbol{e}_2,...$}の定義は以下の通りである。

\boldsymbol{e}_0=    
\begin{pmatrix}
      1 \\
      0 \\
\vdots \\
      0 \\
    \end{pmatrix}
,\boldsymbol{e}_1=    
\begin{pmatrix}
      0 \\
      1 \\
\vdots \\
      0 \\
    \end{pmatrix}
,...,
\boldsymbol{e}_{N-1}=    
\begin{pmatrix}
      0 \\
      0 \\
\vdots \\
      1 \\
    \end{pmatrix}

$\boldsymbol{f}$は別の直交基底ベクトルの組{$\boldsymbol{E}_0,\boldsymbol{E}_1,\boldsymbol{E}_2,...$}と、その係数の組{$F_0,F_1,F_2,...$}で以下のようにも表せる。

\boldsymbol{f} = F_0\boldsymbol{E}_0+F_1\boldsymbol{E}_1+...+F_{N-1}\boldsymbol{E}_{N-1}

ただし、$\boldsymbol{E}_j$ $(j=0,1,...,N-1)$ の定義は以下の通りである。

\boldsymbol{E}_j = 
\begin{pmatrix}
      (e^{i2\pi\frac{0}{N}})^{j} \\
      (e^{i2\pi\frac{1}{N}})^{j}  \\
\vdots \\
      (e^{i2\pi\frac{N-1}{N}})^{j} \\
    \end{pmatrix}

この時、{$\boldsymbol{E}_0,\boldsymbol{E}_1,..$}は以下のような性質を満たす。
$k\neq j$とし、

\boldsymbol{E}_j \cdot \boldsymbol{E}_k^*= N,　
\boldsymbol{E}_j \cdot \boldsymbol{E}_j^*= 0

よって以下のような式が成り立つ。

F_j= \frac{1}{N} \
\boldsymbol{f} \cdot \boldsymbol{E}_j^*
= \frac{1}{N} {}^t\! \ \boldsymbol{E}_j^* \ \boldsymbol{f}

ここで$\boldsymbol{E}_j^T$は$\boldsymbol{E}_j$を転置した横ベクトルである。
またこの係数$F_j$(フーリエ係数)を集めて$\boldsymbol{F}$というベクトルは以下のように定義できる。

\boldsymbol{F}=
    \begin{pmatrix}
      F_0 \\
      F_{1} \\
      \vdots \\
      F_{N-1} \\
    \end{pmatrix}
=
\frac{1}{N} \
\begin{pmatrix}
      {}^t\! \ \boldsymbol{E}_0^* \\
      {}^t\! \ \boldsymbol{E}_1^* \\
      \vdots \\
      {}^t\! \ \boldsymbol{E}_{N-1}^*  \\
    \end{pmatrix} \
\boldsymbol{f}

フーリエ係数のベクトルは元のデータベクトルに上記のような行列をかけることで求められる。

#DFTの行列表現

前述した計算式を簡単に書くために$W_N$という$N \times N$行列を以下のように定義する。

W_N=
\begin{pmatrix}
      {}^t\! \ \boldsymbol{E}_0^* \\
      {}^t\! \ \boldsymbol{E}_1^* \\
      \vdots \\
      {}^t\! \ \boldsymbol{E}_{N-1}^*  \\
    \end{pmatrix} \
=
\frac{1}{N} \
    \begin{pmatrix}
      w_N^{0,0} & w_N^{0,1} & \ldots & w_N^{0,N-1} \\
      w_N^{1,0} & w_N^{1,1} & \ldots & w_N^{0,N-1} \\
      \vdots & \vdots & \ddots & \vdots \\
      w_N^{N-1,0} & w_N^{N-1,1} & \ldots & w_N^{N-1,N-1} \\
    \end{pmatrix}

この時$W_N$の任意の成分$w_N^{j,k}$は以下のように表せる。

w_N^{j,k}=\exp(-i2\pi \frac{j*k}{N}) ,\quad(j,k=0,1,2,...N-1)

この行列を用いてフーリエ係数を表すと以下のようになる。

\boldsymbol{F}=
W_N \
\boldsymbol{f}

#FFTの仕組み
任意のフーリエ係数は元のデータベクトルと、基底ベクトルの複素共役との内積であった。
これを和の形に直して以下の式で表記する。

F_j=\boldsymbol{f} \cdot \boldsymbol{E}_j^* = \frac{1}{N} 
\big\{ \ f_0 w_N^{0,j}+f_1 w_N^{1,j}+..+f_{N-1} w_N^{N-1,j} \ \big\}
=\frac{1}{N} \ \Sigma_{k=0}^{N-1} \ f_k w_N^{k,j}

この和を$k$が奇数の場合と偶数の場合の二通の和の取り方に分ける。

F_j
=\frac{1}{N} \ \Sigma_{k=0}^{N-1} \ f_k w_N^{k,j}
=\frac{1}{N} \big\{ \Sigma_{k'=0}^{ \frac{N}{2} -1} \ f_{2k'} w_N^{2k',j}
+\Sigma_{k'=0}^{ \frac{N}{2} -1} \ f_{2k'+1} w_N^{2k'+1,j}
\big\}
=\frac{1}{N} \big\{ \Sigma_{k'=0}^{ \frac{N}{2} -1} \ f_{2k'} w_N^{2k',j}
+w_N^{1,j} \ \Sigma_{k'=0}^{ \frac{N}{2} -1} \ f_{2k'+1} w_N^{2k',j}
\big\}

ここで

w_{N}^{2k , j} = 
w_{\frac{N}{2}}^{k , j}= 
w_{\frac{N}{2}}^{k , j+\frac{N}{2}}

より

F_{j}
=\frac{1}{N} \big\{ \Sigma_{k'=0}^{ \frac{N}{2} -1} \ f_{k'} w_{\frac{N}{2}}^{k',j}
+w_{N}^{1,j} \ \Sigma_{k'=0}^{ \frac{N}{2} -1} \ f_{k'+1} w_{\frac{N}{2}}^{k',j}
\big\}

上記の式は、元のデータ列から偶数番目の成分のフーリエ変換と奇数番目の成分のフーリエ変換が組み合わさったものである。
ここで、右辺第二項の $w_N^{1 , j} $ 注目すると以下のような性質がある。

w_N^{1 , j+ \frac{N}{2}} = \exp(-i2\pi \frac{j+\frac{N}{2}}{N}) 
= \exp(-i2\pi \frac{j}{N}+-i\pi) 
=- w_N^{1 , j}

つまり以下のような性質がある。

F_{j'}
=\frac{1}{N} \big\{ \Sigma_{k'=0}^{ \frac{N}{2} -1} \ f_{k'} w_{\frac{N}{2}}^{k',j'}
+w_{N}^{1,j'} \ \Sigma_{k'=0}^{ \frac{N}{2} -1} \ f_{k'+1} w_{\frac{N}{2}}^{k',j'}
\big\} 
\\
F_{j'+\frac{N}{2}}
=\frac{1}{N} \big\{ \Sigma_{k'=0}^{ \frac{N}{2} -1} \ f_{k'} w_{\frac{N}{2}}^{k',j'}
-w_{N}^{1,j'} \ \Sigma_{k'=0}^{ \frac{N}{2} -1} \ f_{k'+1} w_{\frac{N}{2}}^{k',j'}
\big\} \\ （j'=0,1,..,\frac{N}{2} -1）

これは大変だ。DFTではバカ正直に$F$の成分一つ一つについて(つまりN個の成分について)計算していたが、FFTでは前半分を計算してしまえば、後半はもうそれを使いまわせばいい。
もちろん掛け算や足し算の計算を余分にしなければならないが、それを考慮しても計算時間が削減できる。
さらに前述したとおり、二つの項はそれぞれ偶数番目、奇数番目の成分のフーリエ変換が組み合わさったものである。このフーリエ変換もFFTでやってしまえばもっともっと計算時間が削減できる。

実際にJuliaとPythonで実装してみてどのようにこのアルゴリズムが動くのか、またどのくらい時間が削減できるかみてみよう。

#Juliaでの実装
まず元データ、実軸、周波数軸の用意から。

using Plots
im=complex(0,1)      #虚数単位 

N=Int(2^6)           #データ数は2の冪乗で用意
X=1                  #xの最大値
dx=X/N               #x軸上のデータ間隔
x=Vector(0:dx:X-dx)  #x軸のデータ

df=1/(X)                      #周波数軸上のデータ間隔
freq  = Vector(0:df:1/dx-df)  #周波数軸の生成
freq[Int(N/2)+1:N] .-= N      #周波数の単位円上での区間を 0~2π から -π~π に

y=sin.(2*pi*3*x)    #波形データを用意

plot(x,
    y,
    xlabel       ="x",   
    ylabel       ="y",            
    fmt=:png,size=(500,300)) #波形をプロット

こんな画像が出力されます。

次にDFTを実装。

function DFT(y)
    l = size(y)[1] 
    F=zeros(Complex,Int(l))
    
    for i in 1: l         #どの周波数を計算するか指定
        for j in 1: l    #シグマをとるためのfor
           F[i] += y[j]*exp(-2*pi*im*((i-1)*(j-1)/l))
        end
    end
    return(F)
end

プロットしてみると以下のようなグラフになる。

F = DFT(y)  #大きさの補正
F ./= (N/2)
F[1] /= 2

plot(freq,
    abs.(F),
    xlabel       ="freq",   
    ylabel       ="amplitude",   
    xlims        =(0,N/2),      
    ylims        =(0,1.1),      
    fmt=:png,size=(500,300))

つぎにFFTを実装してみる。

function W_n(n)         #FFTの第二項で乗算する数値の列を作る。
    X = Vector(0:1:n-1)
    return exp.(-2*pi*im*(X/n))
end

function split_vec(X,n) #入力されたデータを奇数番目の成分と偶数番目の成分に分けるための関数
    arr=[]
    for i in 1:n
        push!(arr,X[i:n:end])
    end
    return arr
end

function FFT(y) 
    N = size(y)[1] 
    F=zeros(Complex,Int(N))
    y_e , y_o=split_vec(y,2) #入力されたデータを奇数番目の成分と偶数番目の成分に分ける
    
    if N < 4  #ある程度まで分割したら、残りはDFTで計算。
        A = DFT(y_e)
        B = W_n(N)[1:Int(N/2)] .* DFT(y_o)
    else
        A = FFT(y_e)
        B = W_n(N)[1:Int(N/2)] .* FFT(y_o)
    end
    
    F[1:Int(N/2)]= A.+B
    F[Int(N/2)+1:end]= A.-B

    return F
end

プロットしてみる。

F = FFT(y)
F ./= (N/2)
F[1] /= 2

plot(freq,
    abs.(F),
    xlabel       ="freq",   
    ylabel       ="amplitude",   
    xlims        =(0,N/2),      
    fmt=:png,size=(500,300))

全く同じグラフが得られた。
一応差をとると、以下のようになる。

plot(freq,
    abs.(FFT(y) .- DFT(y)),
    xlabel       ="freq",   
    ylabel       ="amplitude",   
    xlims        =(0,N/2),  
    ylims        =(0,10^(-12)),  
    fmt=:png,size=(500,300))

これは結果が一致してると言いて良いだろう。

#Juliaによる計算時間の比較。
さて実行時間を比較する。
以下のようなプログラムで計算時間を取得した。

time_DFT=[]  #DFTの計算時間
time_FFT=[]  #FFTの計算時間
n = []　　    #データ数
for i in 2:12
    
    N=Int(2^i)
    push!(n,N)
    
    X=1                  
    dx=X/N              
    x=Vector(0:dx:X-dx)  
    y=sin.(2*pi*3*x)
    push!(time_DFT,@elapsed DFT(y))
    push!(time_FFT,@elapsed FFT(y))
end

以下はFFT、DETの計算時間をプロットしたもの

データ数が多いほどFFTの方が速いことが確かに確認できる。また比をとると以下のようになる。

#Pythonでの実装
まず元データ、実軸、周波数軸の用意から。

import numpy as np
import matplotlib.pyplot as plt
import time

im=complex(0,1)
N=int(2**6)           #データ数は2の冪乗で用意
X=1                  #xの最大値
dx=X/N               #x軸上のデータ間隔
x=np.linspace(0,X-dx,N)  #x軸のデータ
df=1/(X)                      #周波数軸上のデータ間隔
freq  =np.linspace(0,1/dx-df,N)  #周波数軸の生成
freq[int(N*0.5):] -= N

y=np.sin(2*np.pi*3*x)    #波形データを用意

plt.plot(x,y)

plt.title("Graph Title")
plt.xlabel("X")
plt.ylabel("Y")
plt.show()

こんな画像が出力されます。

次にDFTを実装。

def DFT(y):
    N = len(y) 
    F=np.zeros(N,dtype="complex")

    for i in range(N): #どの周波数を計算するか指定
        for j in range(N):  #シグマをとるためのfor
            F[i] += y[j]*np.exp(-2*np.pi*im*((i-1)*(j-1)/N))
    return(F)

プロットしてみると以下のようなグラフになる。

F = DFT(y)  #大きさの補正
F /= (N/2)
F[0] /= 2

plt.plot(freq,np.abs(F))

plt.title("DFT_plot")
plt.xlabel("freq")
plt.ylabel("amplitude")
plt.xlim(0,N/2)
plt.show()

つぎにFFTを実装してみる。

def W_n(n):         #FFTの第二項で乗算する数値の列を作る。
    X = np.linspace(0,N-1,N) 
    return np.exp(-2*np.pi*im*(X/n))

def split_vec(X,n): #入力されたデータを奇数番目の成分と偶数番目の成分に分けるための関数
    arr=[]
    N=len(X)
    for i in range(n):
        arr.append(X[i:N:n])
    return arr

def FFT(y):
    N = len(y)
    F = np.zeros(N,dtype="complex")
    y_e , y_o=split_vec(y,2) #入力されたデータを奇数番目の成分と偶数番目の成分に分ける

    if N < 4 : #ある程度まで分割したら、残りはDFTで計算。
        A = DFT(y_e)
        B = W_n(N)[0:int(N/2)] * DFT(y_o)
    else:
        A = FFT(y_e)
        B = W_n(N)[0:int(N/2)] * FFT(y_o)

    F[0 : int(N/2)] = A + B
    F[int(N/2) : N] = A - B

    return F

プロットしてみる。

F = FFT(y)  #大きさの補正
F /= (N/2)
F[0] /= 2

plt.plot(freq,np.abs(F))

plt.title("Graph Title")
plt.xlabel("X")
plt.ylabel("Y")
plt.xlim(0,N/2)
plt.show()

全く同じグラフが得られた。
これは結果が一致してると言いて良いだろう。

#Pythonによる計算時間の比較。
さて実行時間を比較する。
以下のようなプログラムで計算時間を取得した。

time_DFT=[]  #DFTの計算時間
time_FFT=[]  #FFTの計算時間
n = []              #データ数

for i in range(10):
    N = int(2**(i+2))           #データ数は2の冪乗で用意
    n.append(N)
    X = 1                  #xの最大値
    dx = X/N               #x軸上のデータ間隔
    x = np.linspace(0,X-dx,N)  #x軸のデータ
    y = np.sin(2*np.pi*3*x)    #波形データを用意
    
    t0 = time.time()
    DFT(y)
    time_DFT.append(time.time()-t0)

    t0 = time.time()
    FFT(y)
    time_FFT.append(time.time()-t0)
    
time_DFT = np.array(time_DFT)
time_FFT = np.array(time_FFT)

以下はFFT、DETの計算時間をプロットしたもの

データ数が多いほどFFTの方が速いことが確かに確認できる。また比をとると以下のようになる。

#google colaboratory
それぞれのプログラムは以下のURLのgooglecolaboratoryにノートブック形式で置いてあります。

julia : https://colab.research.google.com/drive/1ng4wHVj1TV5sPvbvRoDOnlpjpDPg0GjS?usp=sharing
python : https://colab.research.google.com/drive/1BSiOckj6xhSJpGvHia76BgAPQ_P5ETCu?usp=sharing

#結論
FFTは、スペクトルを半分計算して残りにその計算結果を流用することで高速化している。
あと、確かにめっちゃ速い。