関数fift2 マニュアル

(The documentation of function fift2)

Last Update: 2021/12/7

◆機能・用途(Purpose)

フーリエスペクトル\(F(f)\)の逆変換 \[\begin{equation} f(t) \equiv \int_{-\infty}^{\infty} F(f) e^{-2\pi ift}df \label{eq.continuous.definition} \end{equation}\] を高速フーリエ変換(FFT)のアルゴリズムを用いて計算する。結果(時系列データ)を長さが2の巾乗の複素数値系列として返す。
Compute the Fourier inverse transformation (\ref{eq.continuous.definition}) of a spectrum \(F(f)\) using the fast Fourier transformation (FFT) algorithm, and return the result (a time series data) as a complex number series of a length which is a power of 2.

◆形式(Format)

#include <sequence/fft.h>
inline struct imsequence2 fift2(struct imsequence2 seq,const double t0)

◆引数(Arguments)

seq

フーリエスペクトル\(F(f)\)を表す構造体。高周波側の半分は低周波側の折り返しの複素共役でなければならない。またメンバsizeは2の巾乗、メンバt0は0.0でなければならない。
A structure which represents the Fourier spectrum \(F(t)\). The upper half frequency components must be the complex conjugates of the lower half with the reversed order. Members size and t0 must be a power of 2 and 0.0, respectively.

逆変換して得られる時系列データの先頭時刻。時系列データをフーリエ変換する際に先頭時刻\(t_0\)の情報が失われるので、逆変換によって完全に元の時系列データに戻すには \(t_0\)の値をスペクトルとは別に与える必要がある。
The beginning time of the transformed time series data. As the beginning time \(t_0\) of a time series data is lost by Fourier transformation, the value of \(t_0\) must be given in addition to the spectrum to completely reconstruct the original time series data by the inverse transformation.

◆戻り値(Return value)

\(F(f)\)をフーリエ逆変換して得られる時系列データ。戻り値のメンバの値は以下のようになる。
A time series data obtained by the Fourier inverse transformation of \(F(f)\). The values of members of the return value are as below.

戻り値のメンバ Member of the return value	値 Value
size	seq.size
t0	引数t0の値 The value of argument t0
dt	(\ref{eq.dt.definition})式の\(\Delta t\) \(\Delta t\) of eq. (\ref{eq.dt.definition})
各\(k\)に対するvalue[k] value[k] for each \(k\)	(\ref{eq.discrete.definition})式の\(f_k\) \(f_k\) 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.

\(N=log_2\)(seq.size)
\(\Delta f=\)seq.dt
\(t_0=\)t0
\(F_n=\)seq.value[n]

◆使用例(Example)

struct sequence f_original;
struct imsequence2 F=fft2_sequence(f);
struct imsequence2 f_reconstructed=fift2(F,f_original.t0);

◆補足(Additional notes)

この関数の戻り値の時系列データは複素数値系列(struct imsequence2型構造体)であり、長さも2の巾乗のままである。戻り値の時系列データを実数値系列(struct sequence型構造体)に直して末尾の不要な部分をカットして必要な長さにする処理まで一括で行う関数として fift2_sequenceを別途用意している。したがってユーザとして逆変換を行うには通常は関数fift2_sequenceを用いれば良い。この関数(fift2)は関数fift2_sequenceの内部で用いられており、ユーザが直接用いるための関数というよりは関数fift2_sequence内部で間接的に用いるための関数という位置づけである。
The return value (time series data) of this function is a complex number series (struct imsequence2-type structure), with the length of a power of 2. There is another, more convenient function fift2_sequence, which converts the time series data to a real number series (struct sequence-type structure) and cut unnecesary parts at the end to make the time series a desired length. Therefore users can usually use the function fift2_sequence. This function (fift2) is used internally within the function fift2_sequence. The main assumed role of this function (fift2) 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.

「Numerical Recipes in C」では時系列データの先頭時刻が0.0であると仮定されている。一方、この関数では時系列データの先頭時刻は任意であり、引数t0で与える。
In “Numerical Recipes in C”, the beginning time of the time series data is assumed to be 0.0. In this function, the beginning time is arbitrary, and is given by an argument t0.
(\ref{eq.continuous.definition})式の積分を台形公式で近似すると和に\(\Delta f\)が掛かった形になる。「Numerical Recipes in C」ではこの\(\Delta f\)が考慮されていない。一方、この関数では和に\(\Delta f\)を掛けることにより、正しい次元を持つ物理量としてのスペクトルから時系列データへの変換を行う。
By approximating the integral in eq. (\ref{eq.continuous.definition}) with the trapezoidal formula, \(\Delta t\) appears as a common factor to the approximated summation form. This \(\Delta t\) is not taken into account in “Numerical Recipes in C”, while it is multiplied with the summation to transform the spectrum to time series data with a correct physical unit in this function.
「Numerical Recipes in C」では時系列データやフーリエスペクトルが1次元実数値配列で表現されるのに対し、この関数では時系列データやフーリエスペクトルを struct imsequence2型構造体を用いて表現する。
The time series data and Fourier spectrum are expressed by 1-D arrays of real numbers in “Numerical Recipes in C”, while they are represented by struct imsequence2-type structures in this function.
「Numerical Recipes in C」ではフーリエスペクトルの配列が時系列データで上書きされる。一方、この関数ではフーリエスペクトルはそのまま残し、新しい構造体を用意して時系列データを返す。
In “Numerical Recipes in C”, the array for the Fourier spectrum is overwritten by the time series data. In this function, the Fourier spectrum is kept unchanged, and the time series data is returned as a new structure.
「Numerical Recipes in C」ではデータサンプル数が2の巾乗になっているか否かのチェックが行われない。一方、この関数ではデータサンプル数をチェックし、 2の巾乗でない場合はプログラムをエラー終了する。
In “Numerical Recipes in C”, the number of data samples is not checked. In this function, a check is conducted whether the number of data samples is a power of 2, and if not, the program finishes as an error.

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

関数fftとの整合が取れるように \[\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}\] \[\begin{equation} \Delta t \equiv \frac{1}{2^N \Delta f} \label{eq.dt.definition} \end{equation}\] と離散化して積分(\ref{eq.continuous.definition})を和で近似すれば \[\begin{eqnarray} f_k &=& \sum_{n=0}^{2^N-1} F_n e^{-2\pi i n\Delta f (t_0+k\Delta t)} \Delta f \nonumber \\ &=& \Delta f \sum_{n=0}^{2^N-1} F_n e^{-2\pi i n \Delta f t_0} e^{-2\pi ink\Delta f \Delta t} \nonumber \\ &=& \Delta f \sum_{n=0}^{2^N-1} F_n e^{-2\pi i n \Delta f t_0} e^{-\frac{2\pi ink}{2^N}} \label{eq.discrete.definition} \end{eqnarray}\] を得る。この関数では、この(\ref{eq.discrete.definition})式で与えられる変換を実行する。
To be consistent with function fft, the time series data and the Fourier spectrum are discretized as eqs. (\ref{eq.discretize.t})-(\ref{eq.dt.definition}). Approximating the integral in eq. (\ref{eq.continuous.definition}) by a summation results in eq. (\ref{eq.discrete.definition}). A conversion given by this equation is computed by this function.

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

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

3-1. \(F_n\)の要素の並べ替え (Sorting the \(F_n\) components)
(\ref{eq.discrete.definition})式の\(n\)についての和の中に登場する \(k\)によらない部分を \[\begin{equation} G_n \equiv F_n e^{-2\pi i n \Delta f t_0} \label{eq.Gn} \end{equation}\] とし、これを\(m\)回並べ替えた配列として \[\begin{equation} G_n^{(m)}=G_{p(m,n)} \hspace{1em} (m=0,1,2,\cdots,N; n=0,1,2,\cdots,2^N-1) \label{eq.Gnm} \end{equation}\] を定義する。ここで\(p(m,n)\)は\(G_n\)の配列要素の並べ替え方を示す配列で、漸化式 \[\begin{equation} p(0,n)=n \label{eq.p0n} \end{equation}\] \[\begin{eqnarray} & & p(m+1,n)= \begin{cases} p(m,n_1(n,l,m)) & (n_{min}^1(l,m) \leq n < n_{max}^1(l,m)) \\ p(m,n_2(n,l,m)) & (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.pmn} \end{eqnarray}\] によって定義されるものとする。ここで \[\begin{equation} n_1(n,l,m)\equiv 2n-l\cdot 2^{N-m} \label{eq.n1} \end{equation}\] \[\begin{equation} n_2(n,l,m)\equiv 2n+1-(l+1)\cdot 2^{N-m} \label{eq.n2} \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}\] とおいた。 \[\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}^2(2^m,m)=2^N \label{eq.nmax2en} \end{equation}\] であるので、\(n=0,1,2,\cdots,2^N-1\)は区間 \([n_{min}^1(0,m),n_{max}^1(0,m)-1]\), \([n_{min}^2(0,m),n_{max}^2(0,m)-1]\), \([n_{min}^1(1,m),n_{max}^1(1,m)-1]\), \([n_{min}^2(1,m),n_{max}^2(1,m)-1]\), \(\cdots\), \([n_{min}^2(2^m,m),n_{max}^2(2^m,m)-1]\)に分割することができる。これらの各区間について(\ref{eq.pmn})式によって\(p(m+1,n)\)を定義すれば、全ての\(n\)に対する\(p(m+1,n)\)を漏れなく重複なく定義される。
Let \(G_n\) be the factor in the summation of eq. (\ref{eq.discrete.definition}) that is independent of \(k\), defined by eq. (\ref{eq.Gn}). Let \(G_n^{(m)}\) be the array \(G_n\) sorted \(m\) times, defined by eq. (\ref{eq.Gnm}), where \(p(m,n)\) is an array that indicates the sorting rule of \(G_n\) defined recursively by eqs. (\ref{eq.p0n}) and (\ref{eq.pmn}). Here, definitions of eqs. (\ref{eq.n1})-(\ref{eq.nmax2}) are used. Because eqs. (\ref{eq.nmax1nmin2})-(\ref{eq.nmax2en}) hold, \(n=0,1,2,\cdots,2^N-1\) can be divided into sections \([n_{min}^1(0,m),n_{max}^1(0,m)-1]\), \([n_{min}^2(0,m),n_{max}^2(0,m)-1]\), \([n_{min}^1(1,m),n_{max}^1(1,m)-1]\), \([n_{min}^2(1,m),n_{max}^2(1,m)-1]\), \(\cdots\), \([n_{min}^2(2^m,m),n_{max}^2(2^m,m)-1]\). Therefore, by defining \(p(m+1,n)\) for each section by eq. (\ref{eq.pmn}), \(p(m+1,n)\) for all \(n\) are defined without missing nor duplication.

3-2. 離散フーリエ変換の計算 (Computations of the discrete Fourier transformation)
\(m=0,1,2,\cdots,N\), \(k=0,1,2,\cdots,2^N-1\)について、 \[\begin{eqnarray} f_k^{(m)} &\equiv& \sum_{n=0}^{2^{N-m}-1} G_{n+l\cdot 2^{N-m}}^{(m)} e^{\frac{-2\pi ink}{2^{N-m}}} \nonumber \\ & & (l\cdot 2^{N-m}\leq k < (l+1)\cdot 2^{N-m}; l=0,1,2\cdots,2^m-1) \label{eq.fkm} \end{eqnarray}\] という量を定義する。\(f_k^{(m)}\)は漸化式 \[\begin{equation} f_k^{(N)}=G_k^{(N)} \label{eq.fkN} \end{equation}\] \[\begin{eqnarray} & & f_k^{(m)}= \begin{cases} f_k^{(m+1)}+E(m)^kf_{k_1(k,m)}^{(m+1)} & (k_{min}^1(l,m) \leq k < k_{max}^1(l,m)) \\ f_{k_2(k,m)}^{(m+1)}+E(m)^kf_k^{(m+1)} & (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.fkm.recursive} \end{eqnarray}\] によって計算できる。ここで \[\begin{equation} k_1(k,m)\equiv k+2^{N-m-1} \label{eq.k1} \end{equation}\] \[\begin{equation} k_2(k,m)\equiv k-2^{N-m-1} \label{eq.k2} \end{equation}\] \[\begin{equation} E(m)\equiv e^{\frac{-2\pi i}{2^{N-m}}} \label{eq.Em} \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}\] とおいた。 \(f_k^{(0)}\)に\(\Delta f\)を掛ければ (\ref{eq.discrete.definition})式の\(f_k\)が得られる。
Let us define \(f_k^{(m)}\) for \(m=0,1,2,\cdots,N\) and \(k=0,1,2,\cdots,2^N-1\) by eq. (\ref{eq.fkm}). The values of \(f_k^{(m)}\) can be computed recursively using eqs. (\ref{eq.fkm}) and (\ref{eq.fkN}), where definitions of eqs. (\ref{eq.k1})-(\ref{eq.kmax2}) are used. Note that \(f_k\) in eq. (\ref{eq.discrete.definition}) is obtained my multiplying \(\Delta f\) to \(f_k^{(0)}\).

3-3. 計算の手順 (Computation steps)
入力スペクトル\(F_n\)を(\ref{eq.Gn})式に代入して\(G_n\)を求める。これを\(G_n^{(0)}\)と見なした上で、漸化式(\ref{eq.p0n})(\ref{eq.pmn})式を繰り返し用いて\(p(N,n)\)を求め、これを(\ref{eq.Gnm})式に代入して\(G_n^{(N)}\)を求める。次に(\ref{eq.fkN})(\ref{eq.fkm.recursive})式を繰り返し用いて \(f_k^{(0)}\)を求め、最後にこれに\(\Delta f\)を掛けて時系列データ\(f_k\)を得る。
The input Fourier spectrum \(F_n\) is inserted into eq. (\ref{eq.Gn}) to compute \(G_n\), which is used as \(G_n^{(0)}\). The values of \(p(N,n)\) are computed by the recursive formula eqs. (\ref{eq.p0n}) and (\ref{eq.pmn}). Inserting \(p(N,n)\) into eq. (\ref{eq.Gnm}) yields \(G_n^{(N)}\). Then eqs. (\ref{eq.fkN}) and (\ref{eq.fkm.recursive}) are used to compute \(f_k^{(0)}\), which is multiplied by \(\Delta f\) to obtain the time series data \(f_k\).