我想使用具有大约 300×200×200
个元素的数组 Intel MKL
来计算3D FFT
. 此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
函数使用的数组的大小? -
如何取出零填充结果并获得实际结果?
这是我对这个问题的尝试,对于任何评论,建议或提示我都会是 absolutely thankful .
#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;
}
1 回答
只有几点提示:
2大小的力量
我不喜欢你计算尺寸的方式
所以让
Nx,Ny,Nz
是输入矩阵的大小和
nx,ny,nz
填充矩阵的大小现在将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中的
是相同的,但您需要添加一个通道和不同的标准化常量
这样你可以有单行缓冲
double lin[2*max(nx,ny,nz)];
并在运行时进行零填充(因此不需要在内存中有更大的矩阵)...
但这涉及应对每个1D DFFT的线路......