5

私は次元 Nx*Ny*Nz の実際の 3D 配列を持っており、 FFTW を使用して各 z 値に対して 2D フーリエ変換を行いたいと考えています。ここで、z インデックスはメモリ内で最も速く変化します。現在、次のコードは期待どおりに機能します。

int Nx = 16; int Ny = 8; int Nz = 3;

// allocate memory
const int dims = Nx * Ny * Nz;

// input data (pre Fourier transform)
double *input = fftw_alloc_real(dims);              

// why is this the required output size?
const int outdims = Nx * (Ny/2 + 1) * Nz;           
// we want to perform the transform out of place
// (so seperate array for output)
fftw_complex *output = fftw_alloc_complex(outdims);    


// setup "plans" for forward and backward transforms
const int rank = 2; const int howmany = Nz;
const int istride = Nz; const int ostride = Nz;
const int idist = 1; const int odist = 1;
int n[] = {Nx, Ny};
int *inembed = NULL, *onembed = NULL;


fftw_plan fp = fftw_plan_many_dft_r2c(rank, n, howmany,
        input, inembed, istride, idist,
        output, onembed, ostride, odist,
        FFTW_PATIENT);

fftw_plan bp = fftw_plan_many_dft_c2r(rank, n, howmany,
        output, onembed, ostride, odist,
        input, inembed, istride, idist,
        FFTW_PATIENT);

私が理解しているように、長さ N の 1D シーケンスを変換するには (N/2 + 1) の複雑な値が必要ですが、代わりにoutdims = (Nx/2 + 1)*(Ny/2 + 1)*Nz2D 変換に期待されるように設定すると、上記のコードがクラッシュするのはなぜですか?

qx = 0 to Nx/2第二に、次を使用して（包括的）からモードの実数部と虚数部にアクセスできると考えているのは正しいです：

#define outputRe(qx,qy,d) ( output[(d) + Nz * ((qy) + (Ny/2 + 1) * (qx))][0] )
#define outputIm(qx,qy,d) ( output[(d) + Nz * ((qy) + (Ny/2 + 1) * (qx))][1] )

EDIT：遊びたい人のための完全なコードとMakefile 。fftw と gsl がインストールされていることを前提としています。

EDIT2：私が正しく理解していれば、インデックス付け（正と負の周波数を許可する）は次のようになります（おそらくマクロには乱雑になりすぎます！）：

#define outputRe(qx,qy,d) ( output[(d) + Nz * ((qy) + (Ny/2 + 1) * ( ((qx) >= 0) ? (qx) : (Nx + (qx)) ) )  ][0] )
#define outputIm(qx,qy,d) ( output[(d) + Nz * ((qy) + (Ny/2 + 1) * ( ((qx) >= 0) ? (qx) : (Nx + (qx)) ) )  ][1] )

for (int qx = -Nx/2; qx < Nx/2; ++qx)
    for (int qy = 0; qy <= Ny/2; ++qy)
        outputRe(qx, qy, d) = ...

whereoutputRe(-Nx/2, qy, d)は と同じデータを指しますoutputRe(Nx/2, qy, d)。実際には、おそらく最初のインデックスをループして周波数に変換するだけで、その逆ではありません!

4

2 に答える 2

5

明確にするために（3Dデータの2D変換に簡単に拡張できるため、2Dに焦点を当てます）：

保管所

Nx * Ny配列にはNx * (Ny / 2 + 1)、フーリエ変換後の複雑な要素が必要です。

まず、y 方向では、負の周波数を複素共役対称性 (実数列の変換から得られる) から再構成できます。y モードは包括的kyに実行されます。したがって、y には複雑な値が必要です。0 to Ny/2Ny/2 + 1

次に、複素数値の y 値に作用しているため、同じ対称性の仮定を使用できない x 方向に変換します。したがって、正と負の周波数を含める必要があるため、x モードは包括的にkx実行されます。-Nx/2 to Nx/2ただしkx = -Nx/2、 とkx = Nx/2は同等であるため、格納されるのは 1 つだけです (こちらを参照)。したがって、x にはNx複雑な値が必要です。

アクセス

tir38 が指摘しているように、変換後の x インデックスは 0 から Nx-1 まで実行されますが、これはモードkxが 0 から Nx-1 まで実行されるという意味ではありません。FFTW は、正の周波数を配列の前半にパックし、次に負の周波数を後半に (逆順で) 次のようにパックします。

kx = 0, 1, 2, ..., Nx/2, -Nx/2 + 1, ..., -2, -1

これらの要素へのアクセスについて考える方法は 2 つあります。最初に tir38 が示唆するように、順番にループしてkx、インデックスからモードを計算できます。

for (int i = 0; i < Nx; i++)
{
    // produces the list of kxs above
    int kx = (i <= Nx/2) ? i : i - Nx;

    // here we index with i, but with the knowledge that the mode is kx
    outputRe(i, ...) = some function of kx 
}

または、モードをループしてkxインデックスに変換できます。

for (int kx = -Nx/2; kx < Nx/2; kx++)
{
    // work out index from mode kx
    int i = (kx >= 0) ? i : Nx + i;

    // here we index with i, but with the knowledge that the mode is kx
    outputRe(i, ...) = some function of kx 
}

2 種類のインデックス作成と残りのコードは、ここにあります。

于 2013-07-03T08:48:31.630 に答える
1

最初の質問:

各 z 値に対して 2D FFT を実行すると、基本的に Nz 個の 2D FFT が実行されます。対称性は 1 次元のみです。単一の [Nx x Ny] 2D FFT の場合、出力は Nx * (Ny/2 +1) になります。したがって、これらの FFT を Nz 回実行しているため、[Nx x (Ny/2 +1+) x Nz] の 3D 出力が得られます。

2番目の質問

はい、それが出力の保存方法です。誰もが Matlab にアクセスできるわけではないことはわかっていますが、私が FFTW を使い始めたとき、私は常に小さな行列を C++ (またはあなたの場合は C) と Matlab の両方と比較していました。

更新 (2 回目の編集に基づく)

それらはまだゼロからインデックス付けされているため、次のようになります。

for (int qx = -Nx/2; qx < Nx/2; ++qx) => for (int gx = 0; gx < Nx; gx++)

対称性は y 軸上にあるため、x データは冗長ではありません。

output(0, a, b) != output(Nx-1, a ,b)

ここで、a、b は y と z の値です。

于 2013-06-28T02:16:13.973 に答える