gpt4 book ai didi

c++ - 使用带零填充的英特尔 MKL 的 3D FFT

转载 作者:塔克拉玛干 更新时间:2023-11-03 07:51:45 24 4
gpt4 key购买 nike

我想使用 Intel MKL 计算一个数组的 3D FFT about 300×200×200 元素。此 3D 数组以列方式存储为 double 类型的一维数组:

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 并防止它被修改(我需要稍后在我的代码中使用它)然后进行反向计算 就地。我也想要零填充。

我的问题是:

  1. 如何执行零填充?
  2. 当计算中包含零填充时,我应该如何处理 FFT 函数使用的数组大小?
  3. 如何取出补零的结果并得到实际结果?

这是我对这个问题的尝试,如果有任何评论、建议或提示,我将非常感谢

#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 填充矩阵的大小

      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 彼此不靠近
  2. 输出很复杂

    • 如果我做对了a是输入实矩阵
    • afft输出复数矩阵
    • 那么为什么不正确分配空间呢?
    • double *afft = new double [2*nx*ny*nz];
    • 复数是实部+虚部,所以每个数有 2 个值
    • 这也适用于结果的最终打印
    • 和一些 "\r\n" 之后的行将有利于查看
  3. 3D 自由傅里叶变换

    • 我既不使用也不了解你们的 DFFT 库
    • 我用的是我自己的,但无论如何 3D DFFT 可以用 1D DFFT 来完成
    • 如果你按行来做...请看这个2D DFCT by 1D DFFT
    • 在 3D 中是相同的,但您需要添加一个 channel 和不同的归一化常数
    • 这样你可以有单行缓冲区double lin[2*max(nx,ny,nz)];
    • 并在运行时进行零填充(因此不需要在内存中有更大的矩阵)...
    • 但这涉及处理每个 1D DFFT 上的线 ...

关于c++ - 使用带零填充的英特尔 MKL 的 3D FFT,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/27329467/

24 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com