Zenn
🕌

NFFT逆変換

2021/12/12に公開
C++
C
#
fft
#
nfft
tech

NFFTとは，

NFFTとは，離散フーリエ変換，

f_j = f(\boldsymbol{x}_j) = \sum_{\boldsymbol{k}} \hat{f}_{\boldsymbol{k}} e^{-2\pi i \boldsymbol{k}\cdot\boldsymbol{x}_j}

に対して，

  • フーリエ級数 \hat{f}_{\boldsymbol{k}} が2N+1個の波数\boldsymbol{k}で与えられている
  • ノード\bm{x}_j \, (j=0,\dots, M-1)，が，不等間隔に与えれている

ときに，元の関数f_jを求めるアルゴリズムのことです．これを行列ベクトルで表すと，

\boldsymbol{f} = \boldsymbol{A} \boldsymbol{\hat{f}},
\boldsymbol{f} = (f_j)_{j=0,\dots, M-1}, \quad \boldsymbol{\hat{f}} = (\hat{f}_{\boldsymbol{k}})_{\boldsymbol{k}}, \quad \boldsymbol{A} = \left(e^{-2\pi i \boldsymbol{k}\cdot\boldsymbol{x}_j } \right)_{j=0,\dots, M-1; \boldsymbol{k}},

となります．

NFFTとそのC言語ライブラリであるNFFTの簡単な説明・使い方は以下の記事を参照してください，

NFFT逆変換

今回の記事では，NFFTの逆変換（INFFT）の説明をします．

NFFTの逆変換（INFFT）とは，ノードとその上の値(\boldsymbol{x}_j, y_j), \, (j=0,\dots,M-1)が与えれた時に，フーリエ級数\boldsymbol{\hat{f}}を求めることです．つまり，以下の線形連立方程式の解\boldsymbol{\hat{f}}を求めることになります．

\boldsymbol{A\hat{f}} = \boldsymbol{y}

等間隔にノード\boldsymbol{x}_jが与れれている通常のFFTとは異なり，フーリエ逆変換が厳密に存在しているとは限りません．フーリエ行列\boldsymbol{A}の随伴行列\boldsymbol{A}^*を左から掛けかけると求まるわけではありません．

ノードの個数がフーリエ級数の個数個数より多い場合と少ない場合で，それぞれ計算法が異なります．

ノードの個数が過剰な場合

ノードの個数が過剰な場合，フーリエ級数\hat{f}が元の関数fを表現できるほど十分ないということになります．

この場合は，最も良いものを最小二乗法で求めます．

実際には，

\min_{\boldsymbol{\hat{f}}} \left(\sum_{j=0}^{M-1} w_j |y_j - f(\boldsymbol{x}_j)|^2\right)^{1/2}, \quad f(\boldsymbol{x}_j) = (\boldsymbol{A\hat{f}})_{j}

を解いているようですが，この重みw_jがどのように決められているのかは，ドキュメントにも書いてありませんでした．判明したら追記するかもしれません．

ノードの個数が不足な場合

ノードの個数が不足している場合は，フーリエ級数\hat{f}が余分に余っている状態です．つまり，適当な\hat{f}が複数個存在する場合があるということです．

この場合は，サンプルyの個数を内挿によって増やしてやる必要があります．
これはすなわち，複数個ある解\hat{f}の中から，最も良さそうなものを選ぶということです．

最も良いものを選ぶ基準は，

\min_{\boldsymbol{\hat{f}}} \left(\sum_{\boldsymbol{k}}\frac{|\hat{f}_{\boldsymbol{k}|}}{\hat{w}_{\boldsymbol{k}}} \right)^{1/2}, \quad\mathrm{subject \,\, to} \quad \boldsymbol{A\hat{f}} = \boldsymbol{y}

と書いてありますが，重み\hat{w}_{\boldsymbol{k}}の詳細は書いてありませんでした．詳細がわかったら追記するかもしれません．

サンプルコード

今回は，コサイン，

\cos(x) = \cos(2\pi x), \quad x \in \left[-\frac{1}{2}, \frac{1}{2}\right]

のフーリエ級数，

\begin{align*} \hat{f}(k) = \begin{cases} 1/2 & \mathrm{for} \quad k = \pm 1 \\ 0 & \mathrm{otherwise} \end{cases} \end{align*}

を求めるコードを示します．

(注意) ドキュメントが3.3でとまっていて，かつ，逆変換のコードは3.4以降でかなり加筆・修正されているようです．ここで示すものも3.5以外では使えなくなっている可能性が高いです．

#include <stdio.h>
#include <math.h>
#include <complex.h>
#include <nfft3.h>

int main(void)
{
    nfft_plan plan;
    solver_plan_complex inverse_plan;
    int N = 4;
    int M = 16;
    int iter = 4;

    nfft_init_1d(&plan, N, M);

    solver_init_complex(&inverse_plan, (nfft_mv_plan_complex *)(&plan));

    nfft_vrand_shifted_unit_double(plan.x, plan.M_total);

    for (int i = 0; i < plan.M_total; i++)
    {
        inverse_plan.y[i] = cos(2 * M_PI * plan.x[i]);
    }

    // initial guess for f_hat
    for (int k = 0; k < plan.N_total; k++)
    {
        inverse_plan.f_hat_iter[k] = 0.0;
    }

    solver_before_loop_complex(&inverse_plan);

    for (int l = 0; l < iter; l++)
    {
        printf("----- %d iteration -----\n", l + 1);

        solver_loop_one_step_complex(&inverse_plan);

        for (int k = 0; k < plan.N_total; k++)
        {
            printf("%e %e\n", creal(inverse_plan.f_hat_iter[k]), cimag(inverse_plan.f_hat_iter[k]));
        }

        printf("\n Residual r=%e\n-----------------------\n", inverse_plan.dot_r_iter);
    }

    solver_finalize_complex(&inverse_plan);
    nfft_finalize(&plan);

    return 0;
}

この結果は以下のようになります，

----- 1 iteration -----
2.967704e-02 1.008778e-01
4.712111e-01 4.425935e-02
4.601135e-03 0.000000e+00
4.712111e-01 -4.425935e-02

 Residual r=2.235322e-01
-----------------------
----- 2 iteration -----
5.654262e-03 8.838023e-03
4.935071e-01 1.039638e-02
-2.747496e-02 -5.693536e-03
5.010603e-01 -1.039638e-02

 Residual r=1.301747e-02
-----------------------
----- 3 iteration -----
3.159839e-04 -1.878513e-03
4.988973e-01 2.774002e-03
1.477890e-04 -4.532437e-04
5.009242e-01 -2.774002e-03

 Residual r=3.316033e-04
-----------------------
----- 4 iteration -----
-1.651681e-17 1.323644e-18
5.000000e-01 -8.873240e-18
-2.119126e-17 -5.392640e-19
5.000000e-01 7.110964e-18

 Residual r=2.474691e-32
-----------------------

コサインだと4回の反復で十分な精度のフーリエ係数が求まります．

以下はコードの詳細な説明です．

初期化

    nfft_plan plan;
    solver_plan_complex inverse_plan;
    int N = 4;
    int M = 16;
    int iter = 4;

    nfft_init_1d(&plan, N, M);

    solver_init_complex(&inverse_plan, (nfft_mv_plan_complex *)(&plan));
  • 内部で，フーリエ変換と擬逆変換(上記説明の逆変換解法)を繰り返すため，
    逆変換用のプランsolver_plan_complexに加え，nfft_planも必要です．
  • 逆変換用プランinverse_planの初期化には，solver_init_complexを使います．引数として，逆変換用プランsolver_planとNFFTプランnfft_planを渡しますが，型としてはnfft_mv_plan_complexに直してから渡す必要があるようです．

値の代入

    nfft_vrand_shifted_unit_double(plan.x, plan.M_total);

    for (int i = 0; i < plan.M_total; i++)
    {
        inverse_plan.y[i] = cos(2 * M_PI * plan.x[i]);
    }

    // initial guess for f_hat
    for (int k = 0; k < plan.N_total; k++)
    {
        inverse_plan.f_hat_iter[k] = 0.0;
    }
  • ここでは，必要な値を代入します．
  • plan.xは，ノードxの値を入れます．今回はランダムに入れています．
  • inverse_plan.yはノードxに対応する関数fの値です．今回はコサインの値を入れています．
  • inverse_plan.f_hat_iterはこれから求める\hat{f}の値の初期値を入れます．反復解法で求めるので，このイニシャルゲスが良ければ，良い値が見つかることになります．

反復前処理

    solver_before_loop_complex(&inverse_plan);
  • solver_before_loop_complexは反復前の処理を実行しています（多分）．詳細はよくわかりません．

反復解法

    for (int l = 0; l < iter; l++)
    {
        solver_loop_one_step_complex(&inverse_plan);
        ...
    }
  • solver_loop_one_step_complexは，実際に反復解法によって，逆変換を求める関数です．１ループごとにこの関数を呼ぶ必要があります．

終了処理

    solver_finalize_complex(&inverse_plan);
    nfft_finalize(&plan);
  • 2つのプランを終了します．

以上です．

GitHubで編集を提案

Discussion