c++ - ゼロパディング付きのインテルMKLを使用した3D FFT

原文 c++ c 3d fft intel-mkl

FFT個の要素を持つ配列のIntel MKLを使用して3D 300×200×200を計算したい。この3D配列は、列ごとにdoubleタイプの1D配列として格納されます。

for( int k = 0; k < nk; k++ ) // Loop through the height.
    for( int j = 0; j < nj; j++ ) // Loop through the rows.
        for( int i = 0; i < ni; i++ ) // Loop through the columns.
        {
            ijk = i + ni * j + ni * nj * k;
            my3Darray[ ijk ] = 1.0;
        }


入力配列でnot-in-place FFTを実行し、それが変更されないようにして(コードで後で使用する必要があります)、逆方向計算in-placeを実行します。また、ゼロパディングも必要です。

私の質問は:


ゼロパディングを実行するにはどうすればよいですか?
計算にゼロパディングが含まれている場合、FFT関数で使用される配列のサイズをどのように処理すればよいですか?
ゼロ詰めの結果を取り出して実際の結果を得るにはどうすればよいですか?


ここに問題への私の試みがあります、私はコメント、提案、またはヒントに絶対に感謝します。

#include <stdio.h>
#include "mkl.h"

int max(int a, int b, int c)
{
     int m = a;
     (m < b) && (m = b); 
     (m < c) && (m = c); 
     return m;
}

void FFT3D_R2C( // Real to Complex 3D FFT.
    double *in, int nRowsIn , int nColsIn , int nHeightsIn ,
    double *out )
{    
    int n  = max( nRowsIn , nColsIn , nHeightsIn  );
    // Round up to the next highest power of 2.
    unsigned int N = (unsigned int) n; // compute the next highest power of 2 of 32-bit n.
    N--;
    N |= N >> 1;
    N |= N >> 2;
    N |= N >> 4;
    N |= N >> 8;
    N |= N >> 16;
    N++;

    /* Strides describe data layout in real and conjugate-even domain. */
    MKL_LONG rs[4], cs[4];

    // DFTI descriptor.
    DFTI_DESCRIPTOR_HANDLE fft_desc = 0;

    // Variables needed for out-of-place computations.
    MKL_Complex16 *in_fft  = new MKL_Complex16 [ N*N*N ];
    MKL_Complex16 *out_fft = new MKL_Complex16 [ N*N*N ];
    double *out_ZeroPadded = new double [ N*N*N ];

    /* Compute strides */
    rs[3] = 1;           cs[3] = 1;
    rs[2] = (N/2+1)*2;   cs[2] = (N/2+1);
    rs[1] = N*(N/2+1)*2; cs[1] = N*(N/2+1);
    rs[0] = 0;           cs[0] = 0;

    // Create DFTI descriptor.
    MKL_LONG sizes[] = { N, N, N };
    DftiCreateDescriptor( &fft_desc, DFTI_DOUBLE, DFTI_REAL, 3, sizes );

    // Configure DFTI descriptor.
    DftiSetValue( fft_desc, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX );
    DftiSetValue( fft_desc, DFTI_PLACEMENT, DFTI_NOT_INPLACE  ); // Out-of-place transformation.
    DftiSetValue( fft_desc, DFTI_INPUT_STRIDES  , rs  );
    DftiSetValue( fft_desc, DFTI_OUTPUT_STRIDES , cs  );

    DftiCommitDescriptor( fft_desc );
    DftiComputeForward  ( fft_desc, in , in_fft  );

    // Change strides to compute backward transform.
    DftiSetValue        ( fft_desc, DFTI_INPUT_STRIDES , cs);
    DftiSetValue        ( fft_desc, DFTI_OUTPUT_STRIDES, rs);
    DftiCommitDescriptor( fft_desc );
    DftiComputeBackward ( fft_desc, out_fft, out_ZeroPadded );

    // Printing the zero padded 3D FFT result.
    for( long long i = 0; i < (long long)N*N*N; i++ )
        printf("%f\n", out_ZeroPadded[i] );

    /* I don't know how to take out the zero padded results and 
       save the actual result in the variable named "out" */

    DftiFreeDescriptor  ( &fft_desc );

    delete[] in_fft;
    delete[] out_ZeroPadded ;
}

int main()
{
    int n = 10;

    double *a    = new double [n*n*n]; // This array is real.
    double *afft = new double [n*n*n]; 

    // Fill the array with some 'real' numbers.
    for( int i = 0; i < n*n*n; i++ )
        a[ i ] = 1.0;

    // Calculate FFT.
    FFT3D_R2C( a, n, n, n, afft );

    printf("FFT results:\n");
    for( int i = 0; i < n*n*n; i++ )
        printf( "%15.8f\n", afft[i] );

    delete[] a;
    delete[] afft;

    return 0;
}
答え
ほんの少しのヒント:


2サイズの威力


サイズの計算方法が気に入らない
したがって、Nx,Ny,Nzを入力行列のサイズとする
およびnx,ny,nzパディングされた行列のサイズ

for (nx=1;nx<Nx;nx<<=1);
for (ny=1;ny<Ny;ny<<=1);
for (nz=1;nz<Nz;nz<<=1);

memsetでゼロパッドを最初にゼロにしてから、行列の行をコピーします
N^3の代わりにnx*ny*nzにパディングすると、大幅に遅くなる可能性があります
nx,ny,nzが互いに近くにない場合

出力が複雑です


うまくいけば、aは入力実数行列です
およびafft出力複素行列
それで、そのためのスペースを正しく割り当てませんか?
double *afft = new double [2*nx*ny*nz];
複素数は実数部と虚数部なので、数値ごとに2つの値
それは結果の最終印刷にも当てはまります
そして、いくつかの"\r\n"の後の行は表示に適しています

3D DFFT


私はあなたのDFFTライブラリを使用したり、知りません
私は自分のものを使用していますが、とにかく3D DFFTは1D DFFTで実行できます
あなたが行でそれをするなら...これを見てください2D DFCT by 1D DFFT
3Dでも同じですが、1つのパスと異なる正規化定数を追加する必要があります
このようにして、単一行バッファdouble lin[2*max(nx,ny,nz)];を持つことができます
実行時にゼロパディングを作成します(メモリに大きな行列を配置する必要はありません)...
しかし、これには各1D DFFTのラインの対処が含まれます...
関連記事

c++ - C++ヘッダーファイルから型のリストを生成する[終了]

c++ - 可変サイズの配列の初期化

c++ - Visual StudioでのLinuxプラットフォーム用のC++コードの記述

c++ - 回転した長方形で画像ROIを抽出する

c++ - プロセス所有者を取得する(Citrix /プロビジョニング)

c++ - SonarQube Visual Studio 2013 C++プラグイン

c++ - Intel GalileoおよびC++ REST SDK

c++ - operator [] const参照の2次元ポインタのオーバーロード

c++ - C++でFLを16進数に変換する

c++ - スレッドC++を使用してシェルコマンドを実行する