LoginSignup
4
5

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

@toya42

趣味でスペクトル法_#002_実FFT

Last updated at Posted at 2019-06-22

概要

　入力あるいは出力が実数のFFT（実FFT）は，データ長が半分の複素FFTに持ち込める．[](さらに，三角変換(Trigonometric Transform)のひとつである離散コサイン変換（discrete cosine transform, DCT）は実FFTを用いて高速に計算が可能である．)文献1を元に，その原理を説明し，Fortranでの実装例を示す．

更新履歴

2019.06.22 : 公開
2019.06.22 : コード部分のインデントを修正

お願い

　本記事は，石岡圭一著『スペクトル法による数値計算入門』1の内容をベースに構成されている．出典よりも詳しい解説を心がけるが，同じアルゴリズムを説明している為どうしても引用に近い形態となる．著作権的に問題があると感じる方がいらっしゃれば，遠慮なく指摘していただきたい．

定義

　実数を下付き$\rm{r}$で表すことにする．また複素数$C$の複素共役は，$C^*$あるいは$\bar{C}$と書く．（#001と$F\leftrightarrow f$など表記に違いがあるため，定義を再掲する）

実DFT（正変換）

　複素DFT（正変換）

\begin{align}
F\left(k\right)=\frac{1}{N}\sum_{x=0}^{N-1} f\left(x\right) \exp\left(-2\pi i \frac{x\cdot k}{N}\right) \tag{1}\\
\left(k=0,1,...,N-1\right)
\end{align}

において，$f(x)=f_{\rm{r}}(x)$．ただし，$N$は偶数とする．

実DFT(逆変換)

　複素DFT（逆変換）

\begin{align}
f\left(x\right)=\sum_{k=0}^{N-1} F\left(k\right) \exp\left(2\pi i \frac{k\cdot x}{N}\right) \tag{2}\\
\left(x=0,1,...,N-1\right)
\end{align}

において，$f(x)=f_{\rm{r}}(x)$．ただし，$N$は偶数とする．

実FFT

　文献1に従って，データ長$N$の実DFTをデータ長$N/2$の複素FFT（逆変換）を利用して計算する方法を述べる．ここで逆変換を用いるのは複素FFTでは，正変換が逆変換のアルゴリズムを用いて計算できるためである．次にDFTの逆変換と正変換の関係式を示す．

\begin{align}
F\left(k\right)
&=\frac{1}{N}\sum_{x=0}^{N-1} f\left(x\right) \exp\left(-2\pi i \frac{x\cdot k}{N}\right) \\
&=\frac{1}{N}\overline{ \sum_{x=0}^{N-1} \overline{f\left(x\right)} \exp\left(2\pi i \frac{x\cdot k}{N}\right)} \tag{3}\\
&\left(k=0,1,...,N-1\right)
\end{align}

これは当然逆も言えて，先に正変換を定義し，逆変換を正変換を用いて表すこともできる（例えばWikipedia2等）．結局，コーディングの際にどちらを選択するかの違いでしかない（が，基本的には先人が作ったライブラリを使用するのでアルゴリズムを学ぶ時以外は気にする必要がない）．

正変換

\begin{align}
F\left(k\right)
=\frac{1}{N}\sum_{x=0}^{N-1} f_{\rm{r}}\left(x\right) \exp\left(-2\pi i \frac{x\cdot k}{N}\right) \tag{4}\\
\left(k=0,1,...,N-1\right)
\end{align}

$f_{\rm{r}}\left(x\right)$が実数なので，$F$について

\begin{align}
F(N-k)
&=\frac{1}{N}\sum_{x=0}^{N-1} f_{\rm{r}}\left(x\right) \exp\left(-2\pi i \frac{x\left(N-k\right)}{N}\right)\\
&=\frac{1}{N}\sum_{x=0}^{N-1} f_{\rm{r}}\left(x\right) \exp\left(2\pi i \frac{x\cdot k}{N}\right)\\
&=F^*\left(k\right)　\tag{5}\\
\end{align}

$\because$ 式(3)より

\begin{align}
F^*\left(k\right)
=\frac{1}{N}\overline{\overline{ \sum_{x=0}^{N-1} \overline{f_{\rm{r}}\left(x\right)} \exp\left(2\pi i \frac{x\cdot k}{N}\right)}}\\
=\frac{1}{N}\sum_{x=0}^{N-1} \overline{f_{\rm{r}}\left(x\right)} \exp\left(2\pi i \frac{x\cdot k}{N}\right) \\
=\frac{1}{N}\sum_{x=0}^{N-1} f_{\rm{r}}\left(x\right) \exp\left(2\pi i \frac{x\cdot k}{N}\right) \tag{6}\\
\end{align}

が成立し，出力の複素数列$F\left(k\right)$で独立なのは，$F(0)$の実部，$F(k),;\left(k=1,2,...,N/2-1\right)$の実部と虚部，$F(N/2)$の実部，の合計$N$個となる．

　長さ$N$の実数入力データ列$f_{\rm{r}}$を使って,長さ$N/2$の複素数データ列$g$に次のように定義する．

\begin{align}
g\left(x\right)=f_{\rm{r}}\left(2x\right)
+ i\; f_{\rm{r}}\left(2x+1\right)\\
\left(x=0,1,...,N/2-1\right)
\end{align}

これを入力として逆変換を考えると，

\begin{align}
G\left(k\right)
&=\sum_{x=0}^{N/2-1}g\left(x\right)\exp\left(2\pi i \frac{x\cdot k}{N/2}\right)\\
&=\sum_{x=0}^{N/2-1}\left[f_{\rm{r}}\left(2x\right)
+ i\; f_{\rm{r}}\left(2x+1\right)\right]\exp\left(2\pi i \frac{x\cdot k}{N/2}\right)\\
&=\sum_{x=0}^{N/2-1}f_{\rm{r}}\left(2x\right)
\exp\left(2\pi i \frac{x\cdot k}{N/2}\right)
+i\sum_{x=0}^{N/2-1}f_{\rm{r}}\left(2x+1\right)\exp\left(2\pi i \frac{x\cdot k}{N/2}\right)\\
&=G_{\rm{even}}(k)+iG_{\rm{odd}}(k)\tag{7}\\
&\left(k=0,...,N/2-1\right)\\
\end{align}

ここで，$G_{\rm{even}},;G_{\rm{odd}}$は次のように定義した．

\begin{align}
G_{\rm{even}}(k)
&=\sum_{x=0}^{N/2-1}f_{\rm{r}}\left(2x\right)
\exp\left(2\pi i \frac{x\cdot k}{N/2}\right)\\
G_{\rm{odd}}(k)
&=\sum_{x=0}^{N/2-1}f_{\rm{r}}\left(2x+1\right)\exp\left(2\pi i \frac{x\cdot k}{N/2}\right)\tag{8}
\end{align}

式(5)の考え方から，$G(N/2-k)$については，

\begin{align}
G\left(N/2-k\right)
&=\sum_{x=0}^{N/2-1}g\left(x\right)\exp\left(2\pi i \frac{x\cdot \left(N/2-k\right)}{N/2}\right)\\
&=\sum_{x=0}^{N/2-1}f_{\rm{r}}\left(2x\right)
\exp\left(-2\pi i \frac{x\cdot k}{N/2}\right)
+i\sum_{x=0}^{N/2-1}f_{\rm{r}}\left(2x+1\right)\exp\left(-2\pi i \frac{x\cdot k}{N/2}\right)\\
&=G^*_{\rm{even}}(k)+iG^*_{\rm{odd}}(k)\tag{9}\\
&\left(k=1,...,N/2-1\right)\\
\end{align}

となる．よって式(7),(9)を連立すれば $1\leq k \leq N/2-1$の範囲で，$G_{\rm{even}},;G_{\rm{odd}}$が得られる．

ところで，求めたい$F$は式(4)で定義されていた．（再掲）

\begin{align}
F\left(k\right)=\frac{1}{N}\sum_{x=0}^{N-1} f_{\rm{r}}\left(x\right) \exp\left(-2\pi i \frac{x\cdot k}{N}\right) \tag{4}\\
\left(k=0,1,...,N-1\right)
\end{align}

これを，$G_{\rm{even}},;G_{\rm{odd}}$を用いて表すことができれば良い．$G_{\rm{even}},;G_{\rm{odd}}$は，それぞれ$f_{\rm{r}}$の偶数項と奇数項を用いて表されていたので．

\begin{align}
F\left(k\right)
&=\frac{1}{N}\sum_{x=0}^{N/2-1}\left[
f_{\rm{r}}\left(2x\right)
\exp\left(-2\pi i \frac{(2x)\cdot k}{N}\right)
+f_{\rm{r}}\left(2x+1\right)
\exp\left(-2\pi i \frac{(2x+1)\cdot k}{N}\right)
\right]\\
&=\frac{1}{N}\sum_{x=0}^{N/2-1}
f_{\rm{r}}\left(2x\right)
\exp\left(-2\pi i \frac{x\cdot k}{N/2}\right)\\
&+\frac{1}{N}\sum_{x=0}^{N/2-1}
f_{\rm{r}}\left(2x+1\right)
\exp\left(-2\pi i \frac{x\cdot k}{N/2}\right)\exp\left(-2\pi i\frac{k}{N}\right)\\
&=\frac{1}{N}\left[G^*_{\rm{even}}(k)+G^*_{\rm{odd}}(k)\exp\left(-2\pi i\frac{k}{N}\right)\right]\tag{10}\\
&\left(k=0,1,...,N-1\right)\\
\end{align}

となる．式(5)より$0\leq k\leq N/2$で$F$が求まれば良い．$G_{\rm{even}},G_{\rm{odd}}$は，$1\leq k\leq N/2$の範囲で式(7),(9)を連立して求めることができる．

\begin{align}
G^*_{\rm{even}}(k)=\frac{1}{2}\left[G(N/2-k)+G^*(k)\right]\\
G^*_{\rm{odd}}(k)=-\frac{i}{2}\left[G(N/2-k)-G^*(k)\right]
\end{align}

よって．$1\leq k\leq N/2$の$F$は，次のように表される．

\begin{align}
F(k)=\frac{1}{2N}
\left[
\left(G(N/2-k)+G^*(k)\right)-i\exp\left(-2\pi i\frac{k}{N}\right)\left(G(N/2-k)-G^*(k)\right)
\right]\tag{11}
\end{align}

最後に，$k=0,N/2$について考える．

  • $k=0$
\begin{align}
F(0) 
=\frac{1}{N}\left[G^*_{\rm{even}}(0)+G^*_{\rm{odd}}(0)\right]\\
=\frac{1}{N}\left[G_{\rm{even}}(0)+G_{\rm{odd}}(0)\right]\\
=\frac{1}{N}\left[\rm{Re}(G(0))+\rm{Im}(G(0))\right]\tag{12}\\
\end{align}

$\because$ $G_{\rm{even}}(0),;G_{\rm{odd}}(0)$は，その定義から実数である．

  • $k=N/2$
\begin{align}
F(N/2) 
&=\frac{1}{N}\sum_{x=0}^{N/2-1}\left[
f_{\rm{r}}\left(2x\right)
\exp\left(-2\pi i \frac{(2x)\cdot N/2}{N}\right)
+f_{\rm{r}}\left(2x+1\right)
\exp\left(-2\pi i \frac{(2x+1)\cdot N/2}{N}\right)
\right]\\
&=\frac{1}{N}\sum_{x=0}^{N/2-1}
f_{\rm{r}}\left(2x\right)
+\frac{1}{N}\sum_{x=0}^{N/2-1}
f_{\rm{r}}\left(2x+1\right)
\exp\left(-2\pi i \frac{x\cdot N/2}{N/2}\right)\exp\left(-2\pi i\frac{N/2}{N}\right)\\
&=\frac{1}{N}\left[G^*_{\rm{even}}(0)-G^*_{\rm{odd}}(0)\right ]\\
&=\frac{1}{N}\left[\rm{Re}(G(0))-\rm{Im}(G(0))\right]\tag{13}\\
\end{align}

よって，$F\left(k\right)$は長さ$N/2$の複素FFTから求められる．

逆変換

\begin{align}
f_{\rm{r}}\left(x\right)=\sum_{k=0}^{N-1} F\left(k\right) \exp\left(2\pi i \frac{k\cdot x}{N}\right) \tag{2}\\
\left(x=0,1,...,N-1\right)
\end{align}

$f_{\rm{r}}$が実数なので，$1\leq k\leq N-1$の$k$について$F(N-k)=F^*(k)$が言える．また$F(N/2)$は実数である．

( $\because$

$F(k)=a_k+ib_k$，$\exp$関数の引数を$\theta_{k,x}$とする．

\begin{align}
f_{\rm{r}}\left(x\right)=\sum_{k=0}^{N-1} \left[(a_k+ib_k)(\cos\theta_{k,x}+i\sin\theta_{k,x})\right]
\end{align}

$\exp$関数の次の関係を踏まえると，

\begin{align}
\exp\left(2\pi i \frac{\left(N-k\right)\cdot x}{N}\right)
=\exp\left(-2\pi i \frac{k\cdot x}{N}\right)\\
=\cos\theta_{k,x}-i\sin\theta_{k,x}
\end{align}

$f_{\rm{r}}$が実数であるための条件は．$F(N-k)=a_k-ib_k=F^*(k)$となる．

また，

\begin{align}
\exp\left(2\pi i \frac{N/2\cdot x}{N}\right)
&=\exp\left(x\pi i \right)\\
&=\cos(x\pi)\\
&=(-1)^x
\end{align}

より，$F(N/2)$は実数である．）

よって，入力複素数列の自由度からこちらも長さ$N/2$の複素FFTに持ち込むことを考える．式(2)を$F(N-k)=F^*(k)$を用いて書き下すと，

\begin{align}
f_{\rm{r}}\left(x\right)
&=\sum_{k=0}^{N-1} F\left(k\right) \exp\left(2\pi i \frac{k\cdot x}{N}\right) \\
&=F(0)
+\sum_{k=0}^{N/2-1}F(k) \exp\left(2\pi i \frac{k\cdot x}{N}\right)\\
&+(-1)^xF(N/2)+\sum_{k=1}^{N/2-1}F(N-k) \exp\left(2\pi i \frac{(N-k)\cdot x}{N}\right)\\
&=F(0)
+\sum_{k=1}^{N/2-1}F(k) \exp\left(2\pi i \frac{k\cdot x}{N}\right)\\
&+(-1)^x\left(F(N/2)+\sum_{k=1}^{N/2-1}F^*(N/2-k) \exp\left(2\pi i \frac{k\cdot x}{N}\right)\right)\tag{14}\\
\end{align}

となる．$(-1)^x$があるので，$x$が偶数と奇数の場合に分けて書くと

\begin{align}
f_{\rm{r}}(2x)
&=F(0)+F(N/2)
+\sum_{k=1}^{N/2-1}F(k) \exp\left(2\pi i \frac{k\cdot x}{N/2}\right)
+\sum_{k=1}^{N/2-1}F^*(N/2-k) \exp\left(2\pi i \frac{k\cdot x}{N/2}\right)\\
f_{\rm{r}}(2x+1)
&=F(0)-F(N/2)
+\exp\left(2\pi i \frac{k}{N}\right)
\left[\sum_{k=1}^{N/2-1}F(k) \exp\left(2\pi i \frac{k\cdot x}{N/2}\right)
-\sum_{k=1}^{N/2-1}F^*(N/2-k) \exp\left(2\pi i \frac{k\cdot x}{N/2}\right)\right]\tag{15}\\
\end{align}

のように表される．$f_{\rm{r}}(2x)+if_{\rm{r}}(2x+1)$を作れば，

\begin{align}
f_{\rm{r}}(2x)+if_{\rm{r}}(2x+1)
&=F(0)+F(N/2)+i\left[F(0)-F(N/2)\right]\\
&+\sum_{k=1}^{N/2-1}
\left[F(k) 
+F^*(N/2-k)+i \exp\left(2\pi i \frac{k}{N}\right)\left(F(k) 
-F^*(N/2-k)\right)\right]
\exp\left(2\pi i \frac{k\cdot x}{N/2}\right)\\
&=\sum_{k=0}^{N/2-1}H(k)\exp\left(2\pi i \frac{k\cdot x}{N/2}\right)\tag{16}\\
\end{align}

となる．結局，複素FFTへの入力複素数列$H$を次のように作れば，

\begin{align}
H(k)=F(k)+F^*(N/2-k)+i \exp\left(2\pi i \frac{k}{N}\right)
\left[F(k) -F^*(N/2-k)\right]\tag{17}\\
\left(k=0,1,...,N/2-1\right)
\end{align}

$f_{\rm{r}}$は出力複素数列の実部と虚部から得られる．

Fortranでの実装例

! Initialization
subroutine fft1d_r_ini(n,itr,w,wf,wb)
    implicit none
    integer(4),intent(in)  :: n			! 実数データ列の長さ
    integer(4),intent(out) :: itr		! log2(n)
    complex(8),dimension(0:n/2-1,n/2),intent(out) :: w			! 回転因子
    complex(8),dimension(0:n/2-1),    intent(out) :: wf,wb
    integer(4) :: k

    call fft1d_c_ini(n/2,itr,w)			! 複素FFTのinitialization subroutine
    do k=0,n/2-1
        wf(k) = iu*exp(-2.0d0*pi*iu*dble(k)/dble(n))
        wb(k) = iu*exp( 2.0d0*pi*iu*dble(k)/dble(n))
    end do

    return
end subroutine fft1d_r_ini

! Forward FFT ( Real to Complex )
subroutine fft1d_r_forward(n,x,itr,a,w,wf)
    implicit none
    integer(4),                 intent(in)  :: n	! データ長
    integer(4),                 intent(in)  :: itr	! log2(n)
    real(8),   dimension(0:n-1),intent(in)  :: x	! 入力実数データ列
    complex(8),dimension(0:n/2),intent(out) :: a	! 出力複素数データ列
    complex(8),dimension(0:n/2-1,n/2),intent(in) :: w	! 回転因子
    complex(8),dimension(0:n/2-1),    intent(in) :: wf
    integer(4) :: j,k
    complex(8),dimension(0:n/2-1) :: y
    real(8) :: ni,c0r,c0i
    complex(8) :: cn2k,cstr

    ni = 1.0d0/dble(n)
    do j=0,n/2-1
        y(j) = x(2*j)+iu*x(2*j+1)				! 複素FFT（逆変換）に入力するデータ列
    end do
    
    call fft1d_c_backward(n/2,y,itr,w)	! 複素FFT（逆変換）
    
    c0r = real(y(0))
    c0i = imag(y(0))
    a(0)   = ni*(c0r+c0i)
    a(n/2) = ni*(c0r-c0i)
    do k=1,n/2-1
        cn2k = y(n/2-k)
	str = conjg(y(k))
	a(k) = 0.50d0*ni*((cn2k+cstr)-wf(k)*(cn2k-cstr))
    end do

    return
end subroutine fft1d_r_forward

! Backward FFT ( Complex to Real )
subroutine fft1d_r_backward(n,a,itr,x,w,wb)
    implicit none
    integer(4),                 intent(in)  :: n,itr
    complex(8),dimension(0:n/2),intent(in)  :: a	! 入力複素数データ列
    real(8),   dimension(0:n-1),intent(out) :: x	! 出力実数データ列
    complex(8),dimension(0:n/2-1,n/2),intent(in) :: w
    complex(8),dimension(0:n/2-1),    intent(in) :: wb
    integer(4) :: j,k
    complex(8),dimension(0:n/2-1) :: d
    complex(8) :: ak,astr

    do k=0,n/2-1
    	ak   = a(k)
    	astr = conjg(a(n/2-k))
        d(k) = (ak+astr)+wb(k)*(ak-astr)
    end do
    
    call fft1d_c_backward(n/2,d,itr,w)	! 複素FFT（逆変換）
    
    do k=0,n/2-1
    	x(2*k)   = real(d(k))
        x(2*k+1) = imag(d(k))
    end do

    return
end subroutine fft1d_r_backward

まとめ

　勉強のため，実FFTのアルゴリズムについてまとめた．入力/出力が実数という制約がある場合のDFTは，出力/入力される独立な複素数はデータ長の半分であり，残りの半分はその複素共役となる．よって，実FFTは正変換も逆変換もデータ長が半分の複素FFTを用いて計算できることを示した．また，Fortranでの実装例も示した．

次回以降の予定

  1. ライブラリ（Intel MKL, ISPACK）の利用方法（+自作コードの性能評価）
  2. 一次元非線形移流方程式の解析（一次元周期境界条件の問題）

  1. 石岡圭一(2004)『スペクトル法による数値計算入門』東京大学出版 ↩ ↩2 ↩3

  2. Wikipedia : 高速フーリエ変換 (https://ja.wikipedia.org/wiki/高速フーリエ変換) ↩

4
5
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
4
5

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?