関数fft2 マニュアル

(The documentation of function fft2)
Last Update: 2021/12/7


◆機能・用途(Purpose)

時系列データ\(f(t)\)のフーリエスペクトル \[\begin{equation} F(f) \equiv \int_{-\infty}^{\infty} f(t) e^{2\pi ift}dt \label{eq.continuous.definition} \end{equation}\] を高速フーリエ変換(FFT)のアルゴリズムを用いて計算する。 時系列データとして長さが2の巾乗で虚部が0.0の複素数値系列 (struct imsequence2型構造体)を与える。
Compute the Fourier spectrum (eq. \ref{eq.continuous.definition}) of a time series data \(f(t)\) using the fast Fourier transformation (FFT) algorithm, where the time series data is given as a complex number series (struct imsequence2-type structure) of a length of a power of 2 with 0.0 values for the imaginary parts.

◆形式(Format)

#include <sequence/fft.h>
inline struct imsequence2 fft2(struct imsequence2 seq)


◆引数(Arguments)

seq 入力時系列データ\(f(t)\)。 虚部が0.0のstruct imsequence2型構造体として与える。 メンバsizeは2の巾乗でなければならない (関数内で長さのチェックは行われる)。
The input time series data \(f(t)\), given as a struct imsequence2-type structure with imaginary parts being 0.0. Member size must be a power of 2; the length is checked within the function.


◆戻り値(Return value)

\(f(t)\)のフーリエスペクトル。戻り値のメンバの値は以下のようになる。
The Fourier spectrum of \(f(t)\). The values of members of the return value are as below.
戻り値のメンバ
Member of the return value 		値
Value
size seq.size
t0 0.0
dt (\ref{eq.df.definition})式の\(\Delta f\)
\(\Delta f\) of eq. (\ref{eq.df.definition})
各\(n\)に対するvalue[n]
value[n] for each \(n\) 		(\ref{eq.discrete.definition})式の\(F_n\)
\(F_n\) of eq. (\ref{eq.discrete.definition})

ここで、(\ref{eq.discrete.definition})式の記号と 引数の対応関係は以下のようになる。
Here, relations between symbols in eq. (\ref{eq.discrete.definition}) and arguments are given as below.


◆使用例(Example)

struct imsequence2 f;
struct imsequence2 F=fft2(f);


◆補足(Additional notes)

この関数では時系列データを長さが2の巾乗の 複素数値系列(struct imsequence2型構造体)として与える必要がある。 \(f(t)\)として任意の長さの実数値の時系列データ(struct sequence型構造体) を直接与えてFFTを行うことができる関数fft2_sequenceを別途用意している。 したがってユーザとしてFFTを行う分には通常は関数fft2_sequenceを用いれば良い。 この関数(fft2)は関数fft2_sequenceの内部で用いられており、 ユーザが直接用いるための関数というよりは 関数fft2_sequence内部で間接的に用いるための関数という位置づけである。
In this function, the time series data must be given as a complex number series (struct imsequence2-type structure) of a length which is a power of 2. There is another, more convenient function fft2_sequence, which can directly use a time series data of real numbers (a struct sequence-type structure) of arbitrary length as \(f(t)\) to perform FFT. Therefore users can usually use the function fft2_sequence for FFT. This function (fft2) is used internally within the function fft2_sequence. The main assumed role of this function (fft2) is this internal call, rather than being directly used by the users.

◆計算式とアルゴリズム (Formula and algorithm)

1. 概要 (Overview)

この関数で用いているアルゴリズムは 「Numerical Recipes in C」(Press et al., 1988, Cambridge University Press) と考え方は同じである。但し、下記の仕様上の違いがある。
The algorithm used in this function is basically same as that used in “Numerical Recipes in C” (Press et al., 1988, Cambridge University Press). However, the following points are different.


2. 離散フーリエ変換 (Discrete Fourier transformation)

時系列データ\(f(t)\)のフーリエ変換は (\ref{eq.continuous.definition})式で与えられる。 \(f(t)\)が有限の時刻範囲において 一定の時間刻み\(\Delta t\)で与えられているものとし、 データサンプル数を\(2^N\)(\(N\):自然数)、 先頭時刻を\(t_0\)とする。 このとき\(f(t)\)の定義域は\(t_0 \leq t < t_0+2^N\Delta t\)となる。 定義域外では\(f(t)=0\)であるものと仮定する。 離散化後の\(f(t)\)および\(F(f)\)をそれぞれ \[\begin{equation} f_k \equiv f(t_0+k\Delta t) \hspace{1em} (k=0,1,2,\cdots,2^N-1) \label{eq.discretize.t} \end{equation}\] \[\begin{equation} F_n \equiv F(n\Delta f) \label{eq.discretize.f} \end{equation}\] と表し、積分(\ref{eq.continuous.definition})を和で近似すれば \[\begin{equation} F_n = \sum_{k=0}^{2^N-1} f_k e^{2\pi i n\Delta f (t_0+k\Delta t)} \Delta t \label{eq.integral.approximated} \end{equation}\] となる。周波数刻みを \[\begin{equation} \Delta f \equiv \frac{1}{2^N \Delta t} \label{eq.df.definition} \end{equation}\] で与えることにすれば \[\begin{equation} F_n = e^{\frac{2\pi i n t_0}{2^N \Delta t}}\Delta t \sum_{k=0}^{2^N-1} f_k e^{\frac{2\pi ink}{2^N}} \hspace{1em} (n=0,1,2,\cdots,2^N-1) \label{eq.discrete.definition} \end{equation}\] を得る。 この関数では(\ref{eq.discrete.definition})式を用いて フーリエ変換を実行する。
The Fourier transformation of a time series data \(f(t)\) is given by eq. (\ref{eq.continuous.definition}). Let \(f(t)\) be given in a finite time range at a constant time interval \(\Delta t\), and let the number of data samples and start time be \(2^N\) (where \(N\) is a natural number) and \(t_0\), respectively. In this case, the definition range of \(f(t)\) is \(t_0\leq t<t_0+2^N\Delta t\). Outside of the definition range, \(f(t)=0\) is assumed. Let the values of \(f(t)\) and \(F(f)\) after discretizations be represented by eqs. (\ref{eq.discretize.t}) and (\ref{eq.discretize.f}), respectively. Then the integration in eq. (\ref{eq.continuous.definition}) can be approximated as (\ref{eq.integral.approximated}). Let the frequency interval be given by eq. (\ref{eq.df.definition}). Then eq. (\ref{eq.discrete.definition}) is obtained. This function calculates the Fourier transformation based on eq. (\ref{eq.discrete.definition}).

3. FFTのアルゴリズム (The algorithm of FFT)

(\ref{eq.discrete.definition})式の計算を高速に行うためのFFTのアルゴリズムは 以下の3ステップから成る。
The FFT algorithm for a fast computation of eq. (\ref{eq.discrete.definition}) is composed of the following three steps.

3-1. \(f_k\)の要素の並べ替え (Sorting the \(f_k\) components)
\(f_k\)の要素を\(m\)回並べ替えた配列として \[\begin{equation} f_k^{(m)}=f_{p(m,k)} \hspace{1em} (m=0,1,2,\cdots,N; k=0,1,2,\cdots,2^N-1) \label{eq.fkm} \end{equation}\] を定義する。ここで\(p(m,k)\)は\(f_k\)の配列要素の並べ替え方を示す配列で、漸化式 \[\begin{equation} p(0,k)=k \label{eq.p0k} \end{equation}\] \[\begin{eqnarray} & & p(m+1,k)= \begin{cases} p(m,k_1(k,l,m)) & (k_{min}^1(l,m) \leq k < k_{max}^1(l,m)) \\ p(m,k_2(k,l,m)) & (k_{min}^2(l,m) \leq k < k_{max}^2(l,m)) \end{cases} \nonumber \\ & & (m\leq N-1; l=0,1,2,\cdots,2^m-1) \label{eq.pmk} \end{eqnarray}\] によって定義されるものとする。ここで \[\begin{equation} k_1(k,l,m)\equiv 2k-l\cdot 2^{N-m} \label{eq.k1} \end{equation}\] \[\begin{equation} k_2(k,l,m)\equiv 2k+1-(l+1)\cdot 2^{N-m} \label{eq.k2} \end{equation}\] \[\begin{equation} k_{min}^1(l,m)\equiv 2l\cdot 2^{N-m-1} \label{eq.kmin1} \end{equation}\] \[\begin{equation} k_{max}^1(l,m)\equiv (2l+1)\cdot 2^{N-m-1} \label{eq.kmax1} \end{equation}\] \[\begin{equation} k_{min}^2(l,m)\equiv (2l+1)\cdot 2^{N-m-1} \label{eq.kmin2} \end{equation}\] \[\begin{equation} k_{max}^2(l,m)\equiv (2l+2)\cdot 2^{N-m-1} \label{eq.kmax2} \end{equation}\] とおいた。 \[\begin{equation} k_{max}^1(l,m)=k_{min}^2(l,m) \label{eq.kmax1kmin2} \end{equation}\] \[\begin{equation} k_{max}^2(l,m)=k_{min}^1(l+1,m) \label{eq.kmax2kmin1} \end{equation}\] \[\begin{equation} k_{min}^1(0,m)=0 \label{eq.kmin1st} \end{equation}\] \[\begin{equation} k_{max}^2(2^m,m)=2^N \label{eq.kmax2en} \end{equation}\] であるので、\(k=0,1,2,\cdots,2^N-1\)は区間 \([k_{min}^1(0,m),k_{max}^1(0,m)-1]\), \([k_{min}^2(0,m),k_{max}^2(0,m)-1]\), \([k_{min}^1(1,m),k_{max}^1(1,m)-1]\), \([k_{min}^2(1,m),k_{max}^2(1,m)-1]\), \(\cdots\), \([k_{min}^2(2^m,m),k_{max}^2(2^m,m)-1]\)に分割することができる。 これらの各区間について(\ref{eq.pmk})式によって\(p(m+1,k)\)を定義すれば、 全ての\(k\)に対する\(p(m+1,k)\)が漏れなく重複なく定義される。
Let \(f_k^{(m)}\) be the array \(f_k\) sorted \(m\) times, defined by eq. (\ref{eq.fkm}), where \(p(m,k)\) is an array that indicates the sorting rule of \(f_k\) defined recursively by eqs. (\ref{eq.p0k}) and (\ref{eq.pmk}). Here, definitions of eqs. (\ref{eq.k1})-(\ref{eq.kmax2}) are used. Because eqs. (\ref{eq.kmax1kmin2})-(\ref{eq.kmax2en}) hold, \(k=0,1,2,\cdots,2^N-1\) can be divided into sections \([k_{min}^1(0,m),k_{max}^1(0,m)-1]\), \([k_{min}^2(0,m),k_{max}^2(0,m)-1]\), \([k_{min}^1(1,m),k_{max}^1(1,m)-1]\), \([k_{min}^2(1,m),k_{max}^2(1,m)-1]\), \(\cdots\), \([k_{min}^2(2^m,m),k_{max}^2(2^m,m)-1]\). Therefore, by defining \(p(m+1,k)\) for each section by eq. (\ref{eq.pmk}), \(p(m+1,k)\) for all \(k\) are defined without missing nor duplication.

3-2. 離散フーリエ変換の計算 (Computations of the discrete Fourier transformation)
\(m=0,1,2,\cdots,N\), \(n=0,1,2,\cdots,2^N-1\)について、 \[\begin{eqnarray} F_n^{(m)} &\equiv& \sum_{k=0}^{2^{N-m}-1} f_{k+l\cdot 2^{N-m}}^{(m)} e^{\frac{2\pi ink}{2^{N-m}}} \nonumber \\ & & (l\cdot 2^{N-m}\leq n < (l+1)\cdot 2^{N-m}; l=0,1,2\cdots,2^m-1) \label{eq.Fnm} \end{eqnarray}\] という量を定義する。\(F_n^{(m)}\)は漸化式 \[\begin{equation} F_n^{(N)}=f_n^{(N)} \label{eq.FnN} \end{equation}\] \[\begin{eqnarray} & & F_n^{(m)}= \begin{cases} F_n^{(m+1)}+E(m)^nF_{n_1(n,m)}^{(m+1)} & (n_{min}^1(l,m) \leq n < n_{max}^1(l,m)) \\ F_{n_2(n,m)}^{(m+1)}+E(m)^nF_n^{(m+1)} & (n_{min}^2(l,m) \leq n < n_{max}^2(l,m)) \end{cases} \nonumber \\ & & (m\leq N-1; l=0,1,2\cdots,2^m-1) \label{eq.Fnm.recursive} \end{eqnarray}\] によって計算できる。ここで \[\begin{equation} n_1(n,m)\equiv n+2^{N-m-1} \label{eq.n1} \end{equation}\] \[\begin{equation} n_2(n,m)\equiv n-2^{N-m-1} \label{eq.n2} \end{equation}\] \[\begin{equation} E(m)\equiv e^{\frac{2\pi i}{2^{N-m}}} \label{eq.Em} \end{equation}\] \[\begin{equation} n_{min}^1(l,m)\equiv 2l\cdot 2^{N-m-1} \label{eq.nmin1} \end{equation}\] \[\begin{equation} n_{max}^1(l,m)\equiv (2l+1)\cdot 2^{N-m-1} \label{eq.nmax1} \end{equation}\] \[\begin{equation} n_{min}^2(l,m)\equiv (2l+1)\cdot 2^{N-m-1} \label{eq.nmin2} \end{equation}\] \[\begin{equation} n_{max}^2(l,m)\equiv (2l+2)\cdot 2^{N-m-1} \label{eq.nmax2} \end{equation}\] とおいた。\(F_n^{(0)}\)は\(f_k\)の離散フーリエ変換になっている。
Let us define \(F_n^{(m)}\) for \(m=0,1,2,\cdots,N\) and \(n=0,1,2,\cdots,2^N-1\) by eq. (\ref{eq.Fnm}). The values of \(F_n^{(m)}\) can be computed recursively using eqs. (\ref{eq.Fnm}) and (\ref{eq.FnN}), where definitions of eqs. (\ref{eq.n1})-(\ref{eq.nmax2}) are used. Note that \(F_n^{(0)}\) is the discrete Fourier transformation of \(f_k\).

3-3. 係数の調整 (Corrections for common factors)
\(F_n^{(0)}\)においては時系列データの先頭時刻が考慮されておらず、 また積分を和で近似した時に登場する\(\Delta t\)も掛かっていない。 すなわち\(F_n^{(0)}\)そのものは正しいフーリエ変換ではない。 フーリエ変換は\(F_n^{(0)}\)を用いて \[\begin{equation} F_n = e^{\frac{2\pi i n t_0}{2^N \Delta t}}\Delta t F_n^{(0)} \label{eq.Fn.factor} \end{equation}\] と計算できる。
In \(F_n^{(0)}\), the beginning time of the time series data and \(\Delta t\) that appears by approximating the integral by summation are not taken into account, meaning that \(F_n^{(0)}\) is not a correct Fourier transformation. The correct Fourier transformation is obtained from \(F_n^{(0)}\) using eq. (\ref{eq.Fn.factor}).

3-4. 計算の手順 (Computation steps)
入力時系列データ\(f_k\)を\(f_k^{(0)}\)と見なした上で、 漸化式(\ref{eq.p0k})(\ref{eq.pmk})式を繰り返し用いて\(p(N,k)\)を求め、 これを(\ref{eq.fkm})式に代入して\(f_k^{(N)}\)を求める。 次に(\ref{eq.FnN})(\ref{eq.Fnm.recursive})式を繰り返し用いて \(F_n^{(0)}\)を求め、最後にこれを(\ref{eq.Fn.factor})式に代入して フーリエ変換\(F_n\)を求める。
The input time series data \(f_k\) is used for \(f_k^{(0)}\). The values of \(p(N,k)\) are computed by recursively applying eqs. (\ref{eq.p0k}) and (\ref{eq.pmk}). Inserting \(p(N,k)\) into eq. (\ref{eq.fkm}) yields \(f_k^{(N)}\). Then eqs. (\ref{eq.FnN}) and (\ref{eq.Fnm.recursive}) are used to compute \(F_n^{(0)}\), which is inserted into eq. (\ref{eq.Fn.factor}) to finally obtain the Fourier transformation \(F_n\).

4. 関数内部でのコーディングに関する説明 (A description for the internal coding in the function)


4-1. 変数名の対応関係 (Relations of the variable names)
関数内における変数名と本マニュアル中の計算式における変数名との対応関係を 下表に示す。
The table below shows the relations between the variable names in the function and those in the formula of this documentation.
関数内における変数名
Variable name in the function 		本マニュアル中の計算式における変数名
Variable name in the formula of this documentation
p[0][k] 偶数の\(m\)に対する\(p(m,k)\)
\(p(m,k)\) for even \(m\)
p[1][k] 奇数の\(m\)に対する\(p(m,k)\)
\(p(m,k)\) for odd \(m\)
∗ptr_p_old_k1 \(p(m,k_1(k,l,m))\)
∗ptr_p_old_k2 \(p(m,k_2(k,l,m))\)
∗ptr_p_new_k \(p(m+1,k)\)
∗ptr_p_new_kmin \(p(m+1,k_{min}^1(l,m))\), \(p(m+1,k_{min}^2(l,m))\)
∗ptr_p_new_kmax \(p(m+1,k_{max}^1(l,m))\), \(p(m+1,k_{max}^2(l,m))\)
F[0][n] 偶数の\(m\)に対する\(F_n^{(m)}\)
\(F_n^{(m)}\) for even \(m\)
F[1][n] 奇数の\(m\)に対する\(F_n^{(m)}\)
\(F_n^{(m)}\) for odd \(m\)
∗ptr_F_old_n \(F_n^{(m+1)}\)
∗ptr_F_old_n1 \(F_{n_1(n,m)}^{(m+1)}\)
∗ptr_F_old_n2 \(F_{n_2(n,m)}^{(m+1)}\)
∗ptr_F_new_n \(F_n^{(m)}\)
∗ptr_F_new_nmin \(n=n_{min}^1(l,m)\), \(n=n_{min}^2(l,m)\)に対する\(F_n^{(m)}\)
\(F_n^{(m)}\) for \(n=n_{min}^1(l,m)\), \(n=n_{min}^2(l,m)\)
∗ptr_F_new_nmax \(n=n_{max}^1(l,m)\), \(n=n_{max}^2(l,m)\)に対する\(F_n^{(m)}\)
\(F_n^{(m)}\) for \(n=n_{max}^1(l,m)\), \(n=n_{max}^2(l,m)\)
powEn[j] \(\exp(2\pi ij/2^N)\)
∗ptr_powEnj \(\exp(2\pi ij/2^N)\)
∗ptr_F0n \(F_n^{(0)}\)
re.value[n] \(F_n\)
∗ptr_re_value_n \(F_n\)


4-2. \(p(N,k)\)の計算コード (A code for computing \(p(N,k)\))
\(p(N,k)\)は(\ref{eq.pmk})-(\ref{eq.kmax2})式を用いて計算できるが、 計算時間増大の原因となる掛け算の多用を避けるために \(k_1(k,l,m)\), \(k_2(k,l,m)\), \(k_{min}^1(l,m)\), \(k_{max}^1(l,m)\), \(k_{min}^2(l,m)\), \(k_{max}^2(l,m)\) の定義式(\ref{eq.k1})-(\ref{eq.kmax2})を直接用いることはせず、 漸化式を用いてこれらの量を計算する。
Although \(p(N,k)\) can be calculated with eqs. (\ref{eq.pmk})-(\ref{eq.kmax2}), direct uses of the definitions of \(k_1(k,l,m)\), \(k_2(k,l,m)\), \(k_{min}^1(l,m)\), \(k_{max}^1(l,m)\), \(k_{min}^2(l,m)\), and \(k_{max}^2(l,m)\), given by eqs. (\ref{eq.k1})-(\ref{eq.kmax2}), increase the computation time because of multiplying operations. To avoid it, recursive formulas are used in the program.
\(k_{min}^1(l,m)\), \(k_{max}^1(l,m)\), \(k_{min}^2(l,m)\), \(k_{max}^2(l,m)\) については漸化式(\ref{eq.kmax1kmin2})-(\ref{eq.kmin1st})および \[\begin{equation} k_{max}^1(l,m)-k_{min}^1(l,m)=k_{max}^2(l,m)-k_{min}^2(l,m)=2^{N-m-1} \label{eq.kmaxmindiff} \end{equation}\] を用いる。 \(k_1(k,l,m)\)については(\ref{eq.k1})(\ref{eq.kmin1})(\ref{eq.kmax1})式より \[\begin{equation} k_1(k_{min}^1(0,m),0,m)=k_1(0,0,m)=0 \label{eq.k1.st} \end{equation}\] \[\begin{equation} k_1(k+1,l,m)=k_1(k,l,m)+2 \label{eq.k1.recursive} \end{equation}\] \[\begin{eqnarray} k_1(k_{min}^1(l+1,m),l+1,m) &=& k_1((2l+2)\cdot 2^{N-m-1},l+1,m) \nonumber \\ &=& 2(2l+2)\cdot 2^{N-m-1}-(l+1)\cdot 2^{N-m} \nonumber \\ &=& (2l+2)\cdot 2^{N-m}-2^{N-m}-l\cdot 2^{N-m} \nonumber \\ &=& (2l+1)\cdot 2^{N-m}-l\cdot 2^{N-m} \nonumber \\ &=& 2(2l+1)\cdot 2^{N-m-1}-l\cdot 2^{N-m} \nonumber \\ &=& k_1((2l+1)\cdot 2^{N-m-1},l,m) \nonumber \\ &=& k_1(k_{max}^1(l,m),l,m) \label{eq.k1.recursive.l} \end{eqnarray}\] が成り立つ。 \(k_2(k,l,m)\)については(\ref{eq.k2})(\ref{eq.kmin2})(\ref{eq.kmax2})式より \[\begin{equation} k_2(k_{min}^2(0,m),0,m)=k_2(2^{N-m-1},0,m)=2\cdot 2^{N-m-1}+1-2^{N-m}=1 \label{eq.k2.st} \end{equation}\] \[\begin{equation} k_2(k+1,l,m)=k_2(k,l,m)+2 \label{eq.k2.recursive} \end{equation}\] \[\begin{eqnarray} k_2(k_{min}^2(l+1,m),l+1,m) &=& k_2((2l+3)\cdot 2^{N-m-1},l+1,m) \nonumber \\ &=& 2(2l+3)\cdot 2^{N-m-1}+1-(l+2)\cdot 2^{N-m} \nonumber \\ &=& (2l+3)\cdot 2^{N-m}+1-2^{N-m}-(l+1)\cdot 2^{N-m} \nonumber \\ &=& (2l+2)\cdot 2^{N-m}+1-(l+1)\cdot 2^{N-m} \nonumber \\ &=& 2(2l+2)\cdot 2^{N-m-1}+1-(l+1)\cdot 2^{N-m} \nonumber \\ &=& k_2((2l+2)\cdot 2^{N-m-1},l,m) \nonumber \\ &=& k_2(k_{max}^2(l,m),l,m) \label{eq.k2.recursive.l} \end{eqnarray}\] が成り立つ。
For \(k_{min}^1(l,m)\), \(k_{max}^1(l,m)\), \(k_{min}^2(l,m)\), and \(k_{max}^2(l,m)\), recursive formulas given by eqs. (\ref{eq.kmax1kmin2})-(\ref{eq.kmin1st}) and (\ref{eq.kmaxmindiff}) are used. For \(k_1(k,l,m)\), eqs. (\ref{eq.k1}), (\ref{eq.kmin1}), and (\ref{eq.kmax1}) yield (\ref{eq.k1.st})-(\ref{eq.k1.recursive.l}). For \(k_1(k,l,m)\), eqs. (\ref{eq.k2}), (\ref{eq.kmin2}), and (\ref{eq.kmax2}) yield (\ref{eq.k2.st})-(\ref{eq.k2.recursive.l}).
これらの漸化式を用いると\(p(N,k)\)は以下のコードで計算することができる。
Using these recursive formulas, \(p(N,k)\) can be computed by the code below.
// コード1(Code 1)
for(m=0;m<N;m++){
     //\(l=0\)のループにおいて最初に用いる\(k_1\), \(k_2\)の設定
     //((\ref{eq.k1.st})(\ref{eq.k2.st})式)
     //Set the initial values of \(k_1\) and \(k_2\) in the loop of \(l=0\)
     //(eqs. (\ref{eq.k1.st}) and (\ref{eq.k2.st}))
     k1=0;
    k2=1;
    kmax=0;
    for(l=0;l<pow(2,m);l++){
         //(\ref{eq.kmax2kmin1})(\ref{eq.kmaxmindiff})式を用いた \(k_{min}^1\), \(k_{max}^1\)の設定と
         //(\ref{eq.pmk})式を用いた \(k_{min}^1 \leq k < k_{max}^1\)に対する\(p(m+1,k)\)の計算、
         //および(\ref{eq.k1.recursive})(\ref{eq.k1.recursive.l})式を用いた \(k_1\)の更新
         //Set \(k_{min}^1\) and \(k_{max}^1\) using eqs. (\ref{eq.kmax2kmin1}) and (\ref{eq.kmaxmindiff}),
         //compute \(p(m+1,k)\) for \(k_{min}^1 \leq k < k_{max}^1\) using eq. (\ref{eq.pmk}),
         //and update \(k_1\) using eqs. (\ref{eq.k1.recursive}) and (\ref{eq.k1.recursive.l})
         kmin=kmax;
        kmax=kmin+pow(2,N-m-1);
        for(k=kmin;k<kmax;k++){
             p[(m+1)%2][k]=p[m%2][k1];
             k1+=2;
        }
         //(\ref{eq.kmax1kmin2})(\ref{eq.kmaxmindiff})式を用いた \(k_{min}^2\), \(k_{max}^2\)の設定と
         //(\ref{eq.pmk})式を用いた \(k_{min}^2 \leq k < k_{max}^2\)に対する\(p(m+1,k)\)の計算、
         //および(\ref{eq.k2.recursive})(\ref{eq.k2.recursive.l})式を用いた \(k_2\)の更新
         //Set \(k_{min}^2\) and \(k_{max}^2\) using eqs. (\ref{eq.kmax1kmin2}) and (\ref{eq.kmaxmindiff}),
         //compute \(p(m+1,k)\) for \(k_{min}^2 \leq k < k_{max}^2\) using eq. (\ref{eq.pmk}),
         //and update \(k_2\) using eqs. (\ref{eq.k2.recursive}) and (\ref{eq.k2.recursive.l})
         kmin=kmax;
        kmax=kmin+pow(2,N-m-1);
        for(k=kmin;k<kmax;k++){
             p[(m+1)%2][k]=p[m%2][k2];
             k2+=2;
        }
    }
}

実際にはp[(m+1)\%2][k]の値を直接変更するのではなくポインタ変数を利用し、 pow(2,m)などを前もって計算して別の変数に格納しておくことにより 計算を高速化している。
In practice, pointer variables are used instead of directly updating p[(m+1)\%2][k] and the exponent values (e.g., pow(2,m)) are computed and stored in variables in advance to make the computation faster.

4-3. \(F_n^{(0)}\)の計算コード (A code for computing \(F_n^{(0)}\))
\(F_n^{(0)}\)は(\ref{eq.Fnm.recursive})-(\ref{eq.nmax2})式を用いて計算できるが、 ここでも掛け算の回避のために(\ref{eq.n1})-(\ref{eq.nmax2})式の直接利用を避け、 漸化式を用いる。
Although \(F_n^{(0)}\) can be calculated with eqs. (\ref{eq.Fnm})-(\ref{eq.nmax2}), direct uses of eqs. (\ref{eq.n1})-(\ref{eq.nmax2}) are avoided again to reduce calculations of products; instead, recursive formulas are used.
\(n_{min}^1\), \(n_{max}^1\), \(n_{min}^2\), \(n_{max}^2\)については \[\begin{equation} n_{max}^1(l,m)=n_{min}^2(l,m) \label{eq.nmax1nmin2} \end{equation}\] \[\begin{equation} n_{max}^2(l,m)=n_{min}^1(l+1,m) \label{eq.nmax2nmin1} \end{equation}\] \[\begin{equation} n_{min}^1(0,m)=0 \label{eq.nmin1st} \end{equation}\] \[\begin{equation} n_{max}^1(l,m)-n_{min}^1(l,m)=n_{max}^2(l,m)-n_{min}^2(l,m)=2^{N-m-1} \label{eq.nmaxmindiff} \end{equation}\] が成り立つ。\(n_1\)については \[\begin{equation} n_1(n_{min}^1(0,m),m)=n_1(0,m)=2^{N-m-1} \label{eq.n1.st} \end{equation}\] \[\begin{equation} n_1(n+1,m)=n_1(n,m)+1 \label{eq.n1.recursive} \end{equation}\] \[\begin{eqnarray} n_1(n_{min}^1(l+1,m),m) &=& n_1((2l+2)\cdot 2^{N-m-1},m) \nonumber \\ &=& (2l+2)\cdot 2^{N-m-1}+2^{N-m-1} \nonumber \\ &=& (2l+1)\cdot 2^{N-m-1}+2^{N-m-1}+2^{N-m-1} \nonumber \\ &=& n_1((2l+1)\cdot 2^{N-m-1},m)+2^{N-m-1} \nonumber \\ &=& n_1(n_{max}^1(l,m),m)+2^{N-m-1} \label{eq.n1.recursive.l} \end{eqnarray}\] が成り立ち、\(n_2\)については \[\begin{equation} n_2(n_{min}^2(0,m),m)=n_2(2^{N-m-1},m)=2^{N-m-1}-2^{N-m-1}=0 \label{eq.n2.st} \end{equation}\] \[\begin{equation} n_2(n+1,m)=n_2(n,m)+1 \label{eq.n2.recursive} \end{equation}\] \[\begin{eqnarray} n_2(n_{min}^2(l+1,m),m) &=& n_2((2l+3)\cdot 2^{N-m-1},m) \nonumber \\ &=& (2l+3)\cdot 2^{N-m-1}-2^{N-m-1} \nonumber \\ &=& (2l+2)\cdot 2^{N-m-1}-2^{N-m-1}+2^{N-m-1} \nonumber \\ &=& n_2((2l+2)\cdot 2^{N-m-1},m)+2^{N-m-1} \nonumber \\ &=& n_2(n_{max}^2(l,m),m)+2^{N-m-1} \label{eq.n2.recursive.l} \end{eqnarray}\] が成り立つ。
For \(n_{min}^1\), \(n_{max}^1\), \(n_{min}^2\), and \(n_{max}^2\), eqs. (\ref{eq.nmax1nmin2})-(\ref{eq.nmaxmindiff}) hold. For \(n_1\) and \(n_2\), eqs. (\ref{eq.n1.st})-(\ref{eq.n1.recursive.l}) and (\ref{eq.n2.st})-(\ref{eq.n2.recursive.l}) hold, respectively.
更に \[\begin{equation} E(m) =e^{\frac{2\pi i}{2^{N-m}}} =e^{\frac{2\pi i\cdot 2^m}{2^N}} \label{eq.Em.arrange} \end{equation}\] \[\begin{equation} E(m)^n =e^{\frac{2\pi in\cdot 2^m}{2^N}} \label{eq.Emn.arrange} \end{equation}\] であるので、\(j=0,1,2,\cdots,2^N-1\)に対する\(e^{\frac{2\pi ij}{2^N}}\)を 前もって配列powEn[j]に代入しておけば \(j_n=n\cdot 2^m\)に対するpowEn[jn]が\(E(m)^n\)を表すことになる。 この\(j_n\)は漸化式 \[\begin{equation} j_0=0, j_{n+1}=j_n+2^m \label{eq.jn.recursive} \end{equation}\] を用いて計算できる。
In addition, eqs. (\ref{eq.Em.arrange}) and (\ref{eq.Emn.arrange}) indicate that if \(e^{\frac{2\pi ij}{2^N}}\) for \(j=0,1,2,\cdots,2^N-1\) are stored in an array powEn[j] in advance, then powEn[jn] for \(j_n=n\cdot 2^m\) represents \(E(m)^n\). The \(j_n\) values can be calculated recursively with eq. (\ref{eq.jn.recursive}).
以上の関係式を用いると \(F_n^{(0)}\)は以下のコードによって計算できることが分かる。
Using these relations, \(F_n^{(0)}\) can be computed by the code below.
// コード2(Code 2)
for(m=N-1;m>=0;m--){
      //\(l=0\)のループにおいて最初に用いる\(n_1\), \(n_2\)の設定
     //((\ref{eq.n1.st})(\ref{eq.n2.st})式)
     //Set the initial values of \(n_1\) and \(n_2\) in the loop of \(l=0\)
     //(eqs. (\ref{eq.n1.st}) and (\ref{eq.n2.st}))
     n1=pow(2,N-m-1);
    n2=0;
    nmax=0;
    for(l=0;l<pow(2,m);l++){
         //(\ref{eq.nmax2nmin1})(\ref{eq.nmaxmindiff})式を用いた \(n_{min}^1\), \(n_{max}^1\)の設定と
         //(\ref{eq.Fnm})式を用いた \(n_{min}^1 \leq n < n_{max}^1\)に対する\(F_n^{(m)}\)の計算、
         //(\ref{eq.jn.recursive})式を用いた\(j_n\)の更新、
         //および(\ref{eq.n1.recursive})(\ref{eq.n1.recursive.l})式を用いた \(n_1\)の更新
         //Set \(n_{min}^1\) and \(n_{max}^1\) using eqs. (\ref{eq.nmax2nmin1}) and (\ref{eq.nmaxmindiff}),
         //compute \(F_n^{(m)}\) for \(n_{min}^1 \leq n < n_{max}^1\) using eq. (\ref{eq.Fnm}),
         //update \(j_n\) using eq. (\ref{eq.jn.recursive}),
         //and update \(n_1\) using eqs. (\ref{eq.n1.recursive}) and (\ref{eq.n1.recursive.l})
         nmin=nmax;
        nmax=nmin+pow(2,N-m-1);
        j=0;
        for(n=nmin;n<nmax;n++){
             F[m%2][n]=F[(m+1)%2][n]+powEn[j]*F[(m+1)%2][n1];
             j+=pow(2,m);
             n1++;
        }
        n1+=pow(2,N-m-1);
         //(\ref{eq.nmax1nmin2})(\ref{eq.nmaxmindiff})式を用いた \(n_{min}^2\), \(n_{max}^2\)の設定と
         //(\ref{eq.Fnm})式を用いた \(n_{min}^2 \leq n < n_{max}^2\)に対する\(F_n^{(m)}\)の計算、
         //(\ref{eq.jn.recursive})式を用いた\(j_n\)の更新、
         //および(\ref{eq.n2.recursive})(\ref{eq.n2.recursive.l})式を用いた \(n_2\)の更新
         //Set \(n_{min}^2\) and \(n_{max}^2\) using eqs. (\ref{eq.nmax1nmin2}) and (\ref{eq.nmaxmindiff}),
         //compute \(F_n^{(m)}\) for \(n_{min}^2 \leq n < n_{max}^2\) using eq. (\ref{eq.Fnm}),
         //update \(j_n\) using eq. (\ref{eq.jn.recursive}),
         //and update \(n_2\) using eqs. (\ref{eq.n2.recursive}) and (\ref{eq.n2.recursive.l})
         nmin=nmax;
        nmax=nmin+pow(2,N-m-1);
        for(n=nmin;n<nmax;n++){
             F[m%2][n]=F[(m+1)%2][n2]+powEn[j]*F[(m+1)%2][n];
             j+=pow(2,m);
             n2++;
        }
        n2+=pow(2,N-m-1);
    }
}

実際にはF[m\%2][n]の値を直接変更するのではなくポインタ変数を利用し、 pow(2,m)などを前もって計算して別の変数に格納しておくことにより 計算を高速化している。
In practice, pointer variables are used instead of directly updating F[m\%2][n] and the exponent values (e.g., pow(2,m)) are computed and stored in variables in advance to make the computation faster.