NFFT逆変換
NFFTとは,
NFFTとは,離散フーリエ変換,
に対して,
- フーリエ級数
が\hat{f}_{\boldsymbol{k}} 個の波数2N+1 で与えられている\boldsymbol{k} - ノード
,が,不等間隔に与えれている\bm{x}_j \, (j=0,\dots, M-1)
ときに,元の関数
となります.
NFFTとそのC言語ライブラリであるNFFTの簡単な説明・使い方は以下の記事を参照してください,
NFFT逆変換
今回の記事では,NFFTの逆変換(INFFT)の説明をします.
NFFTの逆変換(INFFT)とは,ノードとその上の値
等間隔にノード
ノードの個数がフーリエ級数の個数個数より多い場合と少ない場合で,それぞれ計算法が異なります.
ノードの個数が過剰な場合
ノードの個数が過剰な場合,フーリエ級数
この場合は,最も良いものを最小二乗法で求めます.
実際には,
を解いているようですが,この重み
ノードの個数が不足な場合
ノードの個数が不足している場合は,フーリエ級数
この場合は,サンプル
これはすなわち,複数個ある解
最も良いものを選ぶ基準は,
と書いてありますが,重み
サンプルコード
今回は,コサイン,
のフーリエ級数,
を求めるコードを示します.
(注意) ドキュメントが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
は,実際に反復解法によって,逆変換を求める関数です.1ループごとにこの関数を呼ぶ必要があります.
終了処理
solver_finalize_complex(&inverse_plan);
nfft_finalize(&plan);
- 2つのプランを終了します.
以上です.
Discussion