%表題   MATH2 数学処理上位パッケージ
%
%履歴
%\Drireki{
%         90/03/22 塩谷雅人
%         91/12/10 林祥介
%         91/12/20 塩谷雅人
%         92/03/14 塩谷雅人
%         92/04/22 林祥介   (4.1版)
%         94/04/02 塩谷雅人
%        }
%
\Dchapterhead
\label{fftlib}

\section{概要}

  任意の長さのデータについて高速フーリエ変換をおこなう
  サブルーチンパッケージ.
  NCARの数値計算ライブラリより移植した.

  \vspace{1em}
  {\bf サブルーチンの説明}の中の「定義」の節では,
  入出力パラメータの数学的解説に関して次のような表記法をとる:
  処理する配列{\tt X} (長さ{\tt N})の第$i$番目
  $ (i=1, \ldots, {\tt N}) $の
  配列要素について, 入力時の値を$x_{i}$ (小文字),
  出力時の値を$X_{i}$ (大文字)と書く.

  \vspace{1em}
  以下の7つのサブルーチン群の中で初期化をおこなうサブルーチン
  (サブルーチン名が{\tt I}で終わる)は, そのサブルーチン群に属する
  変換ルーチンを用いる際, かならず最初に1回呼ばなければならない.
  ただしそれ以後は, 異なるデータ数を指定するときに限って
  初期化ルーチンを呼べばよい.
  なお, 初期化ルーチンが用いる作業領域は, 同じサブルーチン群に
  属する変換ルーチンを用いている間変更してはならない.
  (この作業領域には, 因数と三角関数表が格納されている).
\section{サブルーチンのリスト}

% \vspace{1em}
  周期実数値データのフーリエ変換をおこなうサブルーチン群.

  \vspace{1em}
  \begin{tabular}{p{5cm}p{10cm}}
    {\tt RFFTI(N,WSAVE)} & 初期化をおこなう. \\
    {\tt RFFTF(N,R,WSAVE)} & フーリエ順変換をおこなう. \\
    {\tt RFFTB(N,R,WSAVE)} & フーリエ逆変換をおこなう.
  \end{tabular}

  \vspace{1em}
  {\tt RFFTI}, {\tt RFFTF}, {\tt RFFTB}の簡易型サブルーチン群.

  \vspace{1em}
  \begin{tabular}{p{5cm}p{10cm}}
    {\tt EZFFTI(N,WSAVE)} & 初期化をおこなう. \\
    {\tt EZFFTF(N,R,A0,A,B,WSAVE)} & フーリエ順変換をおこなう. \\
    {\tt EZFFTB(N,R,A0,A,B,WSAVE)} & フーリエ逆変換をおこなう.
  \end{tabular}

  \vspace{1em}
  奇の周期データのSINE変換をおこなうサブルーチン群.

  \vspace{1em}
  \begin{tabular}{p{5cm}p{10cm}}
    {\tt SINTI(N,WSAVE)} & 初期化をおこなう. \\
    {\tt SINT(N,X,WSAVE)} & SINE変換をおこなう. \\
  \end{tabular}

  \vspace{1em}
  偶の周期データのCOSINE変換をおこなうサブルーチン群.

  \vspace{1em}
  \begin{tabular}{p{5cm}p{10cm}}
    {\tt COSTI(N,WSAVE)} & 初期化をおこなう. \\
    {\tt COST(N,X,WSAVE)} & COSINE変換をおこなう. \\
  \end{tabular}

  \vspace{1em}
  奇数波数成分のみのSIN変換をおこなうサブルーチン群.

  \vspace{1em}
  \begin{tabular}{p{5cm}p{10cm}}
    {\tt SINQI(N,WSAVE)} & 初期化をおこなう. \\
    {\tt SINQF(N,X,WSAVE)} & SINE順変換をおこなう. \\
    {\tt SINQB(N,X,WSAVE)} & SINE逆変換をおこなう.
  \end{tabular}

  \vspace{1em}
  偶数波数成分のみのCOSINE変換をおこなうサブルーチン群.

  \vspace{1em}
  \begin{tabular}{p{5cm}p{10cm}}
    {\tt COSQI(N,WSAVE)} & 初期化をおこなう. \\
    {\tt COSQF(N,X,WSAVE)} & COSINE順変換をおこなう. \\
    {\tt COSQB(N,X,WSAVE)} & COSINE逆変換をおこなう.
  \end{tabular}

  \vspace{1em}
  周期複素数データのフーリエ変換をおこなうサブルーチン群.

  \vspace{1em}
  \begin{tabular}{p{5cm}p{10cm}}
    {\tt CFFTI(N,WSAVE)} & 初期化をおこなう. \\
    {\tt CFFTF(N,C,WSAVE)} & フーリエ順変換をおこなう. \\
    {\tt CFFTB(N,C,WSAVE)} & フーリエ逆変換をおこなう.
  \end{tabular}

\section{サブルーチンの説明}

  \subsection{RFFTI/RFFTF/RFFTB}
  \label{fftlib.sub.rffti}

    \begin{enumerate}
      \item 機能
      \begin{quote}
        周期実数値データのフーリエ変換をおこなう.
        {\tt RFFTI}は初期化をおこなう;
        {\tt RFFTF}はフーリエ順変換をおこなう;
        {\tt RFFTB}はフーリエ逆変換をおこなう.
      \end{quote}
      \item 定義
      \begin{quote}
        $ N $が偶数のとき$ N' = N/2-1 $,
        $ N $が奇数のとき$ N' = (N-1)/2 $とおく. \\
        {\bf 順変換}は次のように定義される.
        \[ R_{1} = \sum_{i=1}^{N}r_{i}, \]
        \[ R_{2k}   =   \sum_{i=1}^{N}r_{i}\cos \frac{2\pi (i-1)k}{N},
           \mbox{\hspace{1em}}
           R_{2k+1} = - \sum_{i=1}^{N}r_{i}\sin \frac{2\pi (i-1)k}{N}
           \mbox{\hspace{1em}}
           ( k = 1, \ldots, N' ). \]
        ただし$ N $が偶数のとき,
        \[ R_{N} = \sum_{i=1}^{N}(-1)^{i-1}r_{i}. \]

        {\bf 逆変換}は次のように定義される. \\
        $ N $が偶数のとき,
        \[ R_{i} = r_{1} + (-1)^{i-1}r_{N}
           + 2 \sum_{k=1}^{N'}
           ( r_{2k}   \cos \frac{2\pi (i-1)k}{N}
           - r_{2k+1} \sin \frac{2\pi (i-1)k}{N})
           \mbox{\hspace{1em}}
           ( i = 1, \ldots, N ). \]
        $ N $が奇数のとき,
        \[ R_{i} = r_{1}
           + 2 \sum_{k=1}^{N'}
           ( r_{2k}   \cos \frac{2\pi (i-1)k}{N}
           - r_{2k+1} \sin \frac{2\pi (i-1)k}{N})
           \mbox{\hspace{1em}}
           ( i = 1, \ldots, N ). \]
      \end{quote}
      \item 呼び出し方法
      \begin{quote}
        {\tt CALL RFFTI(N,WSAVE)}\\
        {\tt CALL RFFTF(N,R,WSAVE)}\\
        {\tt CALL RFFTB(N,R,WSAVE)}
      \end{quote}
      \item パラメーターの説明
      \begin{quote}
        \begin{tabular}{llp{10cm}}
          {\tt N} & {\tt (I)} & 処理するデータの長さ. \\
          {\tt WSAVE} & {\tt (R)} & 作業用配列. 長さは少なくとも
            2{\tt N}+15以上でなければならない. \\
          {\tt R} & {\tt (R)} & 処理する実数型配列.
            入力パラメータでもあり出力パラメータでもある.
            (上記定義参照).
        \end{tabular}
      \end{quote}
      \item 備考
      \begin{enumerate}
        \item この変換では正規化されない.
          つまり{\tt RFFTF}, {\tt RFFTB}を続けて呼ぶと, もとの
          {\tt N}倍の値が返される.
      \end{enumerate}
    \end{enumerate}

  \subsection{EZFFTI/EZFFTF/EZFFTB}
  \label{fftlib.sub.ezffti}

    \begin{enumerate}
      \item 機能
      \begin{quote}
        {\tt RFFTI}, {\tt RFFTF}, {\tt RFFTB}の簡易型サブルーチン.
        {\tt EZFFTI}は初期化をおこなう;
        {\tt EZFFTF}はフーリエ順変換をおこなう;
        {\tt EZFFTB}はフーリエ逆変換をおこなう.
      \end{quote}
      \item 定義
      \begin{quote}
        {\bf 順変換}は次のように定義される.
        ($ N $が偶数のとき$ N' = N/2-1 $,
         $ N $が奇数のとき$ N' = (N-1)/2 $とおく.)
        \[ A_{0} = \frac{1}{N} \sum_{i=1}^{N}r_{i}, \]
        \[ A_{k} = \frac{2}{N}
                   \sum_{i=1}^{N}r_{i}\cos \frac{2\pi (i-1)k}{N},
           \mbox{\hspace{1em}}
           B_{k} = \frac{2}{N}
                   \sum_{i=1}^{N}r_{i}\sin \frac{2\pi (i-1)k}{N}
           \mbox{\hspace{1em}}
           ( k = 1, \ldots, N' ). \]
        ただし$ N $が偶数のとき,
        \[ A_{N/2} = \frac{1}{N}\sum_{i=1}^{N}(-1)^{i-1}r_{i},
           \mbox{\hspace{1em}}
           B_{N/2} = 0. \]

        {\bf 逆変換}は次のように定義される.
        ($ N $が偶数のとき$ N' = N/2 $,
         $ N $が奇数のとき$ N' = (N-1)/2 $とおく.)
        \[ R_{i} = a_{0} + \sum_{k=1}^{N'}
                   ( a_{k} \cos \frac{2\pi (i-1)k}{N}
                   + b_{k} \sin \frac{2\pi (i-1)k}{N})
           \mbox{\hspace{1em}}
           ( i = 1, \ldots, N ). \]
      \end{quote}
      \item 呼び出し方法
      \begin{quote}
        {\tt CALL EZFFTI(N,WSAVE)}\\
        {\tt CALL EZFFTF(N,R,A0,A,B,WSAVE)}\\
        {\tt CALL EZFFTB(N,R,A0,A,B,WSAVE)}
      \end{quote}
      \item パラメーターの説明
      \begin{quote}
        \begin{tabular}{llp{10cm}}
          {\tt N} & {\tt (I)} & 処理するデータの長さ. \\
          {\tt WSAVE} & {\tt (R)} & 作業用配列. 長さは少なくとも
            3{\tt N}+15以上でなければならない. \\
          {\tt R} & {\tt (R)} & 処理する実数型配列.
            {\tt EZFFTF}においては入力パラメータ,
            {\tt EZFFTB}においては出力パラメータである. \\
          {\tt A0} & {\tt (R)} & 上記定義における$ A_{0} $および
            $ a_{0} $. \\
          {\tt A}, {\tt B} & {\tt (R)} & {\tt N}が偶数のとき
            ${\tt N}/2$, {\tt N}が奇数のとき$({\tt N}-1)/2$の長さの
            実数型配列(上記定義参照).
        \end{tabular}
      \end{quote}
      \item 備考
      \begin{enumerate}
        \item この変換では正規化される.
          つまり{\tt EZFFTF}, {\tt EZFFTB}を続けて呼ぶと, もとの
          値が返される.
      \end{enumerate}
    \end{enumerate}

  \subsection{SINTI/SINT}
  \label{fftlib.sub.sinti}

    \begin{enumerate}
      \item 機能
      \begin{quote}
        奇の周期データのSINE変換をおこなう.
        {\tt SINTI}は初期化をおこなう;
        {\tt SINT}はSINE変換をおこなう.
      \end{quote}
      \item 定義
      \begin{quote}
        {\bf 順変換}は次のように定義される.
        \[ X_{k} = 2\sum_{i=1}^{N}x_{i}\sin \frac{\pi ik}{N+1},
           \mbox{\hspace{1em}}
           ( k = 1, \ldots, N ). \]

        {\bf 逆変換}は順変換と同じである.
      \end{quote}
      \item 呼び出し方法
      \begin{quote}
        {\tt CALL SINTI(N,WSAVE)}\\
        {\tt CALL SINT(N,R,WSAVE)}
      \end{quote}
      \item パラメーターの説明
      \begin{quote}
        \begin{tabular}{llp{10cm}}
          {\tt N} & {\tt (I)} & 処理するデータの長さ. \\
          {\tt WSAVE} & {\tt (R)} & 作業用配列. 長さは少なくとも
            2.5{\tt N}+15以上でなければならない. \\
          {\tt R} & {\tt (R)} & 処理する実数型配列.
            入力パラメータでもあり出力パラメータでもある
            (上記定義参照).
        \end{tabular}
      \end{quote}
      \item 備考
      \begin{enumerate}
        \item {\tt SINT}は逆変換でもある.
          また, この変換では正規化されない.
          つまり{\tt SINT}を2回続けて呼ぶと, もとの
          $2({\tt N}+1)$倍の値が返される.
      \end{enumerate}
    \end{enumerate}

  \subsection{COSTI/COST}
  \label{fftlib.sub.costi}

    \begin{enumerate}
      \item 機能
      \begin{quote}
        偶の周期データのCOSINE変換をおこなう.
        {\tt COSTI}は初期化をおこなう;
        {\tt COST}はCOSINE変換をおこなう.
      \end{quote}
      \item 定義
      \begin{quote}
        {\bf 順変換}は次のように定義される.
        \[ X_{k} = x_{1} + (-1)^{k-1}x_{N}
           + 2\sum_{i=2}^{N}x_{i}\cos \frac{\pi (i-1)(k-1)}{N-1},
           \mbox{\hspace{1em}}
           ( k = 1, \ldots, N ). \]

        {\bf 逆変換}は順変換と同じである.
      \end{quote}
      \item 呼び出し方法
      \begin{quote}
        {\tt CALL COSTI(N,WSAVE)}\\
        {\tt CALL COST(N,R,WSAVE)}
      \end{quote}
      \item パラメーターの説明
      \begin{quote}
        \begin{tabular}{llp{10cm}}
          {\tt N} & {\tt (I)} & 処理するデータの長さ. \\
          {\tt WSAVE} & {\tt (R)} & 作業用配列. 長さは少なくとも
            3{\tt N}+15以上でなければならない. \\
          {\tt R} & {\tt (R)} & 処理する実数型配列.
            入力パラメータでもあり出力パラメータでもある
            (上記定義参照).
        \end{tabular}
      \end{quote}
      \item 備考
      \begin{enumerate}
        \item {\tt COST}は逆変換でもある.
          また, この変換では正規化されない.
          つまり{\tt COST}を2回続けて呼ぶと, もとの
          $2({\tt N}-1)$倍の値が返される.
      \end{enumerate}
    \end{enumerate}

  \subsection{SINQI/SINQF/SINQB}
  \label{fftlib.sub.sinqi}

    \begin{enumerate}
      \item 機能
      \begin{quote}
        奇数波数成分のみのSIN変換をおこなう.
       {\tt SINQI}は初期化をおこなう;
       {\tt SINQF}はSINE順変換をおこなう;
       {\tt SINQB}はSINE逆変換をおこなう.
      \end{quote}
      \item 定義
      \begin{quote}
        {\bf 順変換}は次のように定義される.
        \[ X_{k} = (-1)^{k-1}x_{N}
           + 2 \sum_{i=1}^{N-1}r_{i}\sin \frac{\pi (2k-1)i}{2N}
           \mbox{\hspace{1em}}
           ( k = 1, \ldots, N ). \]

        {\bf 逆変換}は次のように定義される.
        \[ X_{i} = 4 \sum_{k=1}^{N}r_{k} \sin \frac{\pi (2k-1)i}{2N}
           \mbox{\hspace{1em}}
           ( i = 1, \ldots, N ). \]
      \end{quote}
      \item 呼び出し方法
      \begin{quote}
        {\tt CALL SINQI(N,WSAVE)}\\
        {\tt CALL SINQF(N,R,WSAVE)}\\
        {\tt CALL SINQB(N,R,WSAVE)}
      \end{quote}
      \item パラメーターの説明
      \begin{quote}
        \begin{tabular}{llp{10cm}}
          {\tt N} & {\tt (I)} & 処理するデータの長さ. \\
          {\tt WSAVE} & {\tt (R)} & 作業用配列. 長さは少なくとも
            3{\tt N}+15以上でなければならない. \\
          {\tt R} & {\tt (R)} & 処理する実数型配列.
            入力パラメータでもあり出力パラメータでもある
            (上記定義参照).
        \end{tabular}
      \end{quote}
      \item 備考
      \begin{enumerate}
        \item この変換では正規化されない.
          つまり{\tt SINQF}, {\tt SINQB}を続けて呼ぶと, もとの
          4{\tt N}倍の値が返される.
      \end{enumerate}
    \end{enumerate}

  \subsection{COSQI/COSQF/COSQB}
  \label{fftlib.sub.cosqi}

    \begin{enumerate}
      \item 機能
      \begin{quote}
        偶数波数成分のみのCOSSINE変換をおこなう.
       {\tt COSQI}は初期化をおこなう;
       {\tt COSQF}はCOSINE順変換をおこなう;
       {\tt COSQB}はCOSINE逆変換をおこなう.
      \end{quote}
      \item 定義
      \begin{quote}
        {\bf 順変換}は次のように定義される.
        \[ X_{k} = x_{1}
           + 2 \sum_{i=2}^{N}r_{i}\cos \frac{\pi (2k-1)(i-1)}{2N}
           \mbox{\hspace{1em}}
           ( k = 1, \ldots, N ). \]

        {\bf 逆変換}は次のように定義される.
        \[ X_{i} = 4 \sum_{k=1}^{N}r_{k} \cos \frac{\pi (2k-1)(i-1)}{2N}
           \mbox{\hspace{1em}}
           ( i = 1, \ldots, N ). \]
      \end{quote}
      \item 呼び出し方法
      \begin{quote}
        {\tt CALL COSQI(N,WSAVE)}\\
        {\tt CALL COSQF(N,R,WSAVE)}\\
        {\tt CALL COSQB(N,R,WSAVE)}
      \end{quote}
      \item パラメーターの説明
      \begin{quote}
        \begin{tabular}{llp{10cm}}
          {\tt N} & {\tt (I)} & 処理するデータの長さ. \\
          {\tt WSAVE} & {\tt (R)} & 作業用配列. 長さは少なくとも
            3{\tt N}+15以上でなければならない. \\
          {\tt R} & {\tt (R)} & 処理する実数型配列.
            入力パラメータでもあり出力パラメータでもある
            (上記定義参照).
        \end{tabular}
      \end{quote}
      \item 備考
      \begin{enumerate}
        \item この変換では正規化されない.
          つまり{\tt COSQF}, {\tt COSQB}を続けて呼ぶと, もとの
          4{\tt N}倍の値が返される.
      \end{enumerate}
    \end{enumerate}

  \subsection{CFFTI/CFFTF/CFFTB}
  \label{fftlib.sub.cffti}

    \begin{enumerate}
      \item 機能
      \begin{quote}
        周期複素数データのフーリエ変換をおこなう.
        {\tt CFFTI}は初期化をおこなう;
        {\tt CFFTF}はフーリエ順変換をおこなう;
        {\tt CFFTB}はフーリエ逆変換をおこなう.
      \end{quote}
      \item 定義
      \begin{quote}
        以下では$ i = \sqrt{-1} $とする.\\
        {\bf 順変換}は次のように定義される.
        \[ C_{k} = \sum_{j=1}^{N}c_{j}
           \exp (-i\frac{2\pi (j-1)(k-1)}{N})
           \mbox{\hspace{1em}}
           ( k = 1, \ldots, N ). \]

        {\bf 逆変換}は次のように定義される.
        \[ C_{j} = \sum_{k=1}^{N}c_{k}
           \exp (i\frac{2\pi (j-1)(k-1)}{N})
           \mbox{\hspace{1em}}
           ( j = 1, \ldots, N ). \]
      \end{quote}
      \item 呼び出し方法
      \begin{quote}
        {\tt CALL CFFTI(N,WSAVE)}\\
        {\tt CALL CFFTF(N,C,WSAVE)}\\
        {\tt CALL CFFTB(N,C,WSAVE)}
      \end{quote}
      \item パラメーターの説明
      \begin{quote}
        \begin{tabular}{llp{10cm}}
          {\tt N} & {\tt (I)} & 処理するデータの長さ. \\
          {\tt WSAVE} & {\tt (R)} & 作業用配列. 長さは少なくとも
            4{\tt N}+15以上でなければならない. \\
          {\tt C} & {\tt (C)} & 処理する複素数型配列.
            入力パラメータでもあり出力パラメータでもある
            (上記定義参照).
        \end{tabular}
      \end{quote}
      \item 備考
      \begin{enumerate}
        \item この変換では正規化されない.
          つまり{\tt CFFTF}, {\tt CFFTB}を続けて呼ぶと, もとの
          {\tt N}倍の値が返される.
      \end{enumerate}
    \end{enumerate}

  \subsection{その他のサブルーチン}
  \label{fftlib.sub.others}

    このパッケージにはこのほかに以下の下位ルーチンがある.
    ここではサブルーチン名をあげるにとどめる.

    \vspace{1em}
    \begin{center}
    \begin{tabular}{llllllll}
      {\tt cfftb1.f} & {\tt cfftf1.f} & {\tt cffti1.f} &
      {\tt cosqb1.f} & {\tt cosqf1.f} & {\tt ezfft1.f} &
      {\tt passb.f } & {\tt passb2.f} \\
      {\tt passb3.f} & {\tt passb4.f} & {\tt passb5.f} &
      {\tt passf.f } & {\tt passf2.f} & {\tt passf3.f} &
      {\tt passf4.f} & {\tt passf5.f} \\
      {\tt pimach.f} & {\tt radb2.f } & {\tt radb3.f } &
      {\tt radb4.f } & {\tt radb5.f } & {\tt radbg.f } &
      {\tt radf2.f } & {\tt radf3.f } \\
      {\tt radf4.f } & {\tt radf5.f } & {\tt radfg.f } &
      {\tt rfftb1.f} & {\tt rfftf1.f} & {\tt rffti1.f} &
      {\tt sint1.f } &
    \end{tabular}
    \end{center}
