matlab - 通过 cuFFT 在逆 FFT 中缩放

Question

每当我绘制使用 cuFFT 的程序获得的值并将结果与​​ Matlab 的结果进行比较时，我都会得到相同形状的图形，并且最大值和最小值的值位于相同的点。然而，cuFFT 产生的值比 Matlab 产生的值大得多。 Matlab代码是

fs = 1000;                              % sample freq
D = [0:1:4]';                           % pulse delay times
t = 0 : 1/fs : 4000/fs;                 % signal evaluation time
w = 0.5;                                % width of each pulse
yp = pulstran(t,D,'rectpuls',w);
filt = conj(fliplr(yp));
xx = fft(yp,1024).*fft(filt,1024);
xx = (abs(ifft(xx)));    

具有相同输入的 CUDA 代码如下:

cufftExecC2C(plan, (cufftComplex *)d_signal, (cufftComplex *)d_signal, CUFFT_FORWARD);
cufftExecC2C(plan, (cufftComplex *)d_filter_signal, (cufftComplex *)d_filter_signal,     CUFFT_FORWARD);
ComplexPointwiseMul<<<blocksPerGrid, threadsPerBlock>>>(d_signal, d_filter_signal, NX);
cufftExecC2C(plan, (cufftComplex *)d_signal, (cufftComplex *)d_signal, CUFFT_INVERSE);

cuFFT 还执行 1024 点 FFT，批量大小为 2。

使用 NX=1024 的比例因子，值不正确。请告诉我该怎么做。

Answer 1

这是将此问题从未回答列表中删除的迟到答案。

您没有提供足够的信息来诊断您的问题，因为您没有指定设置 cuFFT 计划的方式。您甚至没有指定 Matlab 和 cuFFT 的信号是否具有完全相同的形状(因此您只有缩放比例)或者您具有大约相同的形状。但是，让我提出以下两个意见:

1. yp 向量有 4000 个元素；与 fft(yp,1024) 相反，您通过将信号截断为 1024 元素来执行 FFT；
2. 逆 cuFFT 不按向量元素的数量执行缩放。

为了方便起见(它可能对其他用户有用)，我在下面报告了一个简单的 FFT-IFFT 方案，其中还包括使用 CUDA Thrust 库执行的缩放。

#include <cufft.h>
#include <thrust/host_vector.h>
#include <thrust/device_vector.h>

/*********************/
/* SCALE BY CONSTANT */
/*********************/
class Scale_by_constant
{
    private:
        float c_;

    public:
        Scale_by_constant(float c) { c_ = c; };

        __host__ __device__ float2 operator()(float2 &a) const
        {
            float2 output;

            output.x = a.x / c_;
            output.y = a.y / c_;

            return output;
        }

};

int main(void){

    const int N=4;

    // --- Setting up input device vector    
    thrust::device_vector<float2> d_vec(N,make_cuComplex(1.f,2.f));

    cufftHandle plan;
    cufftPlan1d(&plan, N, CUFFT_C2C, 1);

    // --- Perform in-place direct Fourier transform
    cufftExecC2C(plan, thrust::raw_pointer_cast(d_vec.data()),thrust::raw_pointer_cast(d_vec.data()), CUFFT_FORWARD);

    // --- Perform in-place inverse Fourier transform
    cufftExecC2C(plan, thrust::raw_pointer_cast(d_vec.data()),thrust::raw_pointer_cast(d_vec.data()), CUFFT_INVERSE);

    thrust::transform(d_vec.begin(), d_vec.end(), d_vec.begin(), Scale_by_constant((float)(N)));

    // --- Setting up output host vector    
    thrust::host_vector<float2> h_vec(d_vec);

    for (int i=0; i<N; i++) printf("Element #%i; Real part = %f; Imaginary part: %f\n",i,h_vec[i].x,h_vec[i].y);

    getchar();
}