c++ - FFTW 与 Matlab FFT

Question

我在 matlab central 上发布了这个，但没有得到任何回复，所以我想我会在这里重新发布。

我最近在 Matlab 中编写了一个在 for 循环中使用 FFT 的简单例程； FFT 在计算中占主导地位。出于实验目的，我在 mex 中编写了相同的例程，它调用了 FFTW 3.3 库。事实证明，对于非常大的数组，matlab 例程比 mex 例程运行得更快(大约快两倍)。 mex 例程使用智慧并执行相同的 FFT 计算。我也知道 matlab 使用 FFTW，但他们的版本是否可能稍微优化一些？我什至使用了 FFTW_EXHAUSTIVE 标志，对于大型数组，它的速度仍然比 MATLAB 对应的慢两倍。此外，我确保我使用的 matlab 是带有“-singleCompThread”标志的单线程，并且我使用的 mex 文件未处于 Debug模式。只是好奇是否是这种情况-或者是否有一些我不知道的 matlab 正在使用的优化。谢谢。

这是墨西哥的部分:

void class_cg_toeplitz::analysis() {
// This method computes CG iterations using FFTs
    // Check for wisdom
    if(fftw_import_wisdom_from_filename("cd.wis") == 0) {
        mexPrintf("wisdom not loaded.\n");
    } else {
        mexPrintf("wisdom loaded.\n");
    }

    // Set FFTW Plan - use interleaved FFTW
    fftw_plan plan_forward_d_buffer;    
    fftw_plan plan_forward_A_vec;       
    fftw_plan plan_backward_Ad_buffer;
    fftw_complex *A_vec_fft;
    fftw_complex *d_buffer_fft;
    A_vec_fft = fftw_alloc_complex(n);
    d_buffer_fft = fftw_alloc_complex(n);

    // CREATE MASTER PLAN - Do this on an empty vector as creating a plane 
    // with FFTW_MEASURE will erase the contents; 
    // Use d_buffer
    // This is somewhat dangerous because Ad_buffer is a vector; but it does not
    // get resized so &Ad_buffer[0] should work
    plan_forward_d_buffer = fftw_plan_dft_r2c_1d(d_buffer.size(),&d_buffer[0],d_buffer_fft,FFTW_EXHAUSTIVE);
    plan_forward_A_vec = fftw_plan_dft_r2c_1d(A_vec.height,A_vec.value,A_vec_fft,FFTW_WISDOM_ONLY);
    // A_vec_fft.*d_buffer_fft will overwrite d_buffer_fft
    plan_backward_Ad_buffer = fftw_plan_dft_c2r_1d(Ad_buffer.size(),d_buffer_fft,&Ad_buffer[0],FFTW_EXHAUSTIVE);

    // Get A_vec_fft
    fftw_execute(plan_forward_A_vec);

    // Find initial direction - this is the initial residual
    for (int i=0;i<n;i++) {
        d_buffer[i] = b.value[i];
        r_buffer[i] = b.value[i];
    }    

    // Start CG iterations
    norm_ro = norm(r_buffer);
    double fft_reduction = (double)Ad_buffer.size(); // Must divide by size of vector because inverse FFT does not do this
    while (norm(r_buffer)/norm_ro > relativeresidual_cutoff) {        
        // Find Ad - use fft
        fftw_execute(plan_forward_d_buffer);    
        // Get A_vec_fft.*fft(d) - A_vec_fft is only real, but d_buffer_fft
        // has complex elements; Overwrite d_buffer_fft        
        for (int i=0;i<n;i++) {
            d_buffer_fft[i][0] = d_buffer_fft[i][0]*A_vec_fft[i][0]/fft_reduction;
            d_buffer_fft[i][1] = d_buffer_fft[i][1]*A_vec_fft[i][0]/fft_reduction;
        }        
        fftw_execute(plan_backward_Ad_buffer); 

        // Calculate r'*r
        rtr_buffer = 0;
        for (int i=0;i<n;i++) {
            rtr_buffer = rtr_buffer + r_buffer[i]*r_buffer[i];
        }    

        // Calculate alpha
        alpha = 0;
        for (int i=0;i<n;i++) {
            alpha = alpha + d_buffer[i]*Ad_buffer[i];
        }    
        alpha = rtr_buffer/alpha;

        // Calculate new x
        for (int i=0;i<n;i++) {
            x[i] = x[i] + alpha*d_buffer[i];
        }   

        // Calculate new residual
        for (int i=0;i<n;i++) {
            r_buffer[i] = r_buffer[i] - alpha*Ad_buffer[i];
        }   

        // Calculate beta
        beta = 0;
        for (int i=0;i<n;i++) {
            beta = beta + r_buffer[i]*r_buffer[i];
        }  
        beta = beta/rtr_buffer;

        // Calculate new direction vector
        for (int i=0;i<n;i++) {
            d_buffer[i] = r_buffer[i] + beta*d_buffer[i];
        }  

        *total_counter = *total_counter+1;
        if(*total_counter >= iteration_cutoff) {
            // Set total_counter to -1, this indicates failure
            *total_counter = -1;
            break;
        }
    }

    // Store Wisdom
    fftw_export_wisdom_to_filename("cd.wis");

    // Free fft alloc'd memory and plans
    fftw_destroy_plan(plan_forward_d_buffer);
    fftw_destroy_plan(plan_forward_A_vec);
    fftw_destroy_plan(plan_backward_Ad_buffer);
    fftw_free(A_vec_fft);
    fftw_free(d_buffer_fft);
};

这是 matlab 部分:

% Take FFT of A_vec.
A_vec_fft = fft(A_vec); % Take fft once

% Find initial direction - this is the initial residual 
x = zeros(n,1); % search direction
r = zeros(n,1); % residual
d = zeros(n+(n-2),1); % search direction; pad to allow FFT
for i = 1:n
    d(i) = b(i); 
    r(i) = b(i); 
end

% Enter CG iterations
total_counter = 0;
rtr_buffer = 0;
alpha = 0;
beta = 0;
Ad_buffer = zeros(n+(n-2),1); % This holds the product of A*d - calculate this once per iteration and using FFT; only 1:n is used
norm_ro = norm(r);

while(norm(r)/norm_ro > 10^-6)
    % Find Ad - use fft
    Ad_buffer = ifft(A_vec_fft.*fft(d)); 

    % Calculate rtr_buffer
    rtr_buffer = r'*r;

    % Calculate alpha    
    alpha = rtr_buffer/(d(1:n)'*Ad_buffer(1:n));

    % Calculate new x
    x = x + alpha*d(1:n);

    % Calculate new residual
    r = r - alpha*Ad_buffer(1:n);

    % Calculate beta
    beta = r'*r/(rtr_buffer);

    % Calculate new direction vector
    d(1:n) = r + beta*d(1:n);      

    % Update counter
    total_counter = total_counter+1; 
end

就时间而言，对于 N = 50000 和 b = 1:n，mex 大约需要 10.5 秒，matlab 大约需要 4.4 秒。我正在使用 R2011b。谢谢

Answer 1

由于我不了解 MATLAB FFT 实现的任何细节，因此仅提供一些观察而不是明确的答案:

• 根据您的代码，我可以看到速度差异的两种解释:
  • 速度差异可以通过 FFT 优化级别的差异来解释
  • MATLAB 中的 while 循环执行次数要少得多

我假设您已经研究了第二个问题并且迭代次数是可比较的。 (如果不是，这很可能是一些准确性问题，值得进一步调查。)

现在，关于 FFT 速度比较:

• 是的，理论上 FFTW 比其他高级 FFT 实现更快，但只有在将苹果与苹果进行比较时才有意义:在这里，您正在比较更底层的实现，在汇编级别，其中不仅算法的选择，而且针对特定处理器和具有不同技能的软件开发人员的实际优化都在起作用
• 一年来，我已经在许多处理器上优化或审查了优化的 FFT 组装(我从事基准测试行业)，出色的算法只是故事的一部分。有一些非常特定于您正在编码的架构的考虑因素(考虑延迟、指令调度、优化寄存器使用、内存中的数据排列、考虑分支采用/未采用延迟等)并且会产生差异与算法的选择同样重要。
• 在 N=500000 时，我们还讨论了大内存缓冲区:这是另一扇门，可以进行更多优化，可以很快地针对您运行代码的平台进行调整:您在避免缓存未命中方面做得如何算法决定了数据的流动方式以及软件开发人员可能使用了哪些优化来有效地将数据导入和导出内存。
• 虽然我不知道 MATLAB FFT 实现的细节，但我很确定 DSP 工程师大军已经(并且仍在)优化它，因为它是众多设计的关键。这很可能意味着 MATLAB 拥有合适的开发人员组合来生成更快的 FFT。