関数fft2_pow5_analytical マニュアル

(The documentation of function fft2_pow5_analytical)
Last Update: 2024/4/24


◆機能・用途(Purpose)

時間関数pow5 のフーリエスペクトルを解析式を用いて計算する。
Compute the Fourier spectrum of a time function pow5 using an analytical formula.
時間関数pow5は\(\tau_p\)を正の定数として \[\begin{equation} f(t)= \begin{cases} 0 & (t<0) \\ 10\left(\frac{t}{\tau_p}\right)^3 -15\left(\frac{t}{\tau_p}\right)^4 +6\left(\frac{t}{\tau_p}\right)^5 & (0\leq t \leq \tau_p) \\ 1 & (\tau_p<t) \end{cases} \label{eq.timefunc} \end{equation}\] で定義され、そのフーリエスペクトルは \[\begin{eqnarray} F_{-}(\omega) &\equiv & \int_{-\infty}^{\infty}f(t)e^{-i\omega t}dt \nonumber \\ &=& \frac{60\tau_p}{(-i\omega\tau_p)^4}\left(1-e^{-i\omega\tau_p}\right) +\frac{360\tau_p}{(-i\omega\tau_p)^5}\left(1+e^{-i\omega\tau_p}\right) +\frac{720\tau_p}{(-i\omega\tau_p)^6}\left(1-e^{-i\omega\tau_p}\right) \label{eq.spectrum.negative} \end{eqnarray}\] となる ( 関数pow5のマニュアル参照)。 但しこれはフーリエ変換の定義に\(e^{-i\omega t}\)を用いた場合の式であり、 関数fft2に合わせて\(e^{+i\omega t}\)を用いる場合には (\ref{eq.spectrum.negative})式の\(\omega\)を\(-\omega\)で置き換えて \[\begin{eqnarray} F_{+}(\omega) &\equiv & \int_{-\infty}^{\infty}f(t)e^{+i\omega t}dt \nonumber \\ &=& \frac{60\tau_p}{(i\omega\tau_p)^4}\left(1-e^{i\omega\tau_p}\right) +\frac{360\tau_p}{(i\omega\tau_p)^5}\left(1+e^{i\omega\tau_p}\right) +\frac{720\tau_p}{(i\omega\tau_p)^6}\left(1-e^{i\omega\tau_p}\right) \label{eq.spectrum.positive} \end{eqnarray}\] となる。この関数では(\ref{eq.spectrum.positive})式を用いて フーリエスペクトルを計算する。
A time function pow5 is defined by Eq. (\ref{eq.timefunc}), where \(\tau_p\) is a positive constant. The Fourier spectrum of it is given as Eq. (\ref{eq.spectrum.negative}) (see the documentation of function create_timefunc_pow5) if \(e^{-i\omega t}\) is used for the definition of the spectrum. If \(e^{+i\omega t}\) is used instead (to be consistent with function fft2), the formula for the corresponding Fourier spectrum is obtained by replacing \(\omega\) in Eq. (\ref{eq.spectrum.negative}) with \(-\omega\), and the result is Eq. (\ref{eq.spectrum.positive}). This equation (\ref{eq.spectrum.positive}) is used in this function to compute the Fourier spectrum.
但し、(\ref{eq.spectrum.positive})式の分母に\(\omega\)が登場するので \(\omega=0\)の成分(DC成分)は計算できない。 そこでDC成分については 関数fft2 と等価な計算式を用いて数値的に計算する。関数fft2で用いている計算式は \[\begin{equation} F_n \equiv F_{-}(n\Delta f) =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}} \label{eq.fft2} \end{equation}\] である (\(\Delta f\):周波数刻み、 \(\Delta t\):時間刻み、 \(t_0\):時系列データの先頭時刻、 \(2^N\):時系列データのサンプル数、 \(f_k\):\(k\)番目のサンプル時刻における時系列データの値)。 \(\omega=0\)の成分は(\ref{eq.fft2})式で\(n=0\)として \[\begin{equation} F_0=\Delta t\sum_{k=0}^{2^N-1}f_k \label{eq.DC} \end{equation}\] と計算できる。
However, the Fourier spectrum at \(\omega=0\) (the DC component) cannot be computed analytically because of \(\omega\) in the denominator of Eq. (\ref{eq.spectrum.positive}). Therefore, the DC component is computed numerically using an equation equivalent to that used in function fft2 that is Eq. (\ref{eq.fft2}), where \(\Delta f\) is a frequency step, \(\Delta t\) is a time step, \(t_0\) is the beginning time of the input time series data, \(2^N\) is the number of samples of the time series data, and \(f_k\) is the value of the time series data at \(k\)th time sample. The component for \(\omega=0\) can be computed by letting \(n=0\) in Eq. (\ref{eq.fft2}) and the result is Eq. (\ref{eq.DC}).

◆形式(Format)

#include <sequence/fft.h>
inline struct imsequence2 fft2_pow5_analytical
(const struct sequence timefunc,const double tp)


◆引数(Arguments)

timefunc フーリエ変換したい時間関数。 関数create_timefunc を用いて作成したpow5関数でなければならない。
A time function to convert to the Fourier spectrum. It must be a pow5 function created by function create_timefunc.
tp 時間関数の\(\tau_p\)の値。 引数timefuncで与える時間関数に用いたものと同じでなければならない。
The value of \(\tau_p\) of the time function, which must be same as that used for argument timefunc.


◆戻り値(Return value)

引数timefuncで与えた時間関数のフーリエスペクトル。 (\ref{eq.spectrum.positive})(\ref{eq.DC})式を用いて計算される。 ナイキスト周波数から高周波側は 関数fill_spectrum2_upper_half を用いて計算される。
The Fourier spectrum of the time function given by argument timefunc, computed by Eqs. (\ref{eq.spectrum.positive}) and (\ref{eq.DC}). The spectral components above the Nyquist frequency are computed by function fill_spectrum2_upper_half.

◆使用例(Example)

struct sequence timefunc=create_timefunc ("pow5",10.0,0.0,100000,0.0,0.01,0);
struct imsequence2 spectrum=fft2_pow5_analytical (timefunc,10.0);

赤字で示した2つの\(\tau_p\)が共通でなければならない。
The two \(\tau_p\) shown by red must be equal.

◆検証(Validation)

数値積分で得られるフーリエスペクトルとの比較により検証した。 この検証には0.01秒刻み、長さ1000秒、時定数\(\tau_p=10\) sのpow5関数を用いた。 そのフーリエスペクトルをこの関数(解析式)を用いて計算した場合と 関数fft2_sequence (数値積分)を用いて計算した場合とで同じ結果が得られるかをチェックした。
The output of this function was examined by comparison with a Fourier spectrum obtained by numerical integration. This examination used a pow5 function of 1000 s long with a time constant \(\tau_p=10\) s sampled at every 0.01 s. Its Fourier spectra computed by this function (based on an analytical equation) and by 関数fft2_sequence (numerical integral) were compared to check the results.
数値積分においてはpow5関数そのものを用いると正しい結果が得られない。 これはpow5関数が\(t\rightarrow\infty\)で0に収束しないことによる。 すなわちpow5関数のフーリエ変換は \[\begin{equation} F(\omega) =\int_0^{\tau_p} \left[10\left(\frac{t}{\tau_p}\right)^3 -15\left(\frac{t}{\tau_p}\right)^4 +6\left(\frac{t}{\tau_p}\right)^5\right] e^{i\omega t}dt +\int_{\tau_p}^{\infty}e^{i\omega t}dt \label{eq.spectrum.pow5} \end{equation}\] となり、この第2項の積分を有限の時刻で打ち切ることによって無視できない誤差が生じる。 そこでpow5関数の導関数のフーリエ変換 \[\begin{equation} F_d(\omega) =\int_0^{\tau_p} \left[\frac{30\left(\frac{t}{\tau_p}\right)^2 -60\left(\frac{t}{\tau_p}\right)^3 +30\left(\frac{t}{\tau_p}\right)^4}{\tau_p}\right] e^{i\omega t}dt \label{eq.spectrum.pow5_diff} \end{equation}\] を数値的に計算し、次いで\(F(\omega)=F_d(\omega)/(-i\omega)\)の関係を用いて \(F(\omega)\)を計算する。実際の手順は
  1. 関数create_timefuncを用いてpow5関数の導関数を作成する
  2. そのフーリエスペクトルを関数fft2_sequenceを用いて計算する
  3. 得られたフーリエスペクトルの各周波数成分を\(-i\omega\)で割る
という3ステップである。
The numerical integral of a pow5 function does not yield a correct Fourier spectrum. This is because the pow5 function does not converge to zero for \(t\rightarrow\infty\): the Fourier transformation of the pow5 function is given by Eq. (\ref{eq.spectrum.pow5}), and truncating the integral of the 2nd term of this equation at a finite time results in a significant error. To avoid this problem, the Fourier transformation of the derivative of a pow5 function (F_d(\omega); Eq. \ref{eq.spectrum.pow5_diff}) was computed, and then a relation \(F(\omega)=F_d(\omega)/(-i\omega)\) was applied to obtain the Fourier spectrum \(F(\omega)\). This procedure is composed of the following three steps:
  1. create the derivative of a pow5 function using function create_timefunc;
  2. compute its Fourier spectrum using function fft2_sequence; and
  3. divide each spectral component of the result by \(-i\omega\).
一方、解析式に基づくフーリエスペクトルは
  1. 関数create_timefuncを用いてpow5関数を作成する
  2. そのフーリエスペクトルを関数fft2_pow5_analyticalを用いて計算する
という2ステップで得られる。
The Fourier spectrum based on an analytical equation is computed by the following two steps:
  1. create a pow5 function using function create_timefunc; and
  2. compute its Fourier spectrum using function fft2_pow5_analytical.
下図はこのようにして計算したフーリエスペクトルを示している。 青緑が数値積分(関数fft2_sequence使用)、 赤線が解析式(関数fft2_pow5_analytical使用)である。 両者がよく一致していることからこの関数が正しい結果を与えることが分かる。
The figure below shows the Fourier spectra computed by these two different approaches; the blue-green line is the result of a numerical integration calculated with a function fft2_sequence, and the red line shows an analytical result calculated with a function fft2_pow5_analytical. The excellent fit between the two lines indicates the validity of this function.