使用cuFFT求解傅里叶变换时的最优尺寸选择

前言

傅里叶变换（Fourier Transform）是数字信号处理领域中一个非常重要的数学变换方法，用来实现信号从时域到频域的变换过程。

离散傅里叶变换（Discrete Fourier Transform，DFT）是连续傅里叶变换在离散系统中的表示形式，由于其计算量大，在很长一段时间内其应用受到了极大限制。

快速傅里叶变换（Fast Fourier Transform，FFT），是快速计算序列的离散傅里叶变换（DFT）或其逆变换的方法。是由Cooley和Tukey于20世纪60年代（1965年）合作发表之后开始为人所知。它能够将计算DFT的复杂度从只用DFT定义计算需要的 $O(n^2)$降低到$O(n\log n)$，其中 $n$ 为数据大小。大大提高了DFT的运算速度，从而使DFT在实际应用中得到了快速发展。

对于FFT的CPU实现，一般可以通过调用FFTW库实现。

而对于FFT的GPU实现，一般通过调用cuFFT库实现。

最优尺寸选择

在cuFFT的文档中明确说道：

注意第一条：

Algorithms highly optimized for input sizes that can be written in the
form $2^a\times 3^b\times 5^c\times 7^d$. In general the smaller the prime factor, the better the performance, i.e., powers of two are fastest.

也就是说要使用cuFFT达到比较高效的计算性能，待计算数据的尺寸应该为质数2、3、5、7的幂次方乘积。当尺寸仅为质数2的幂次方时，计算速度最快。（注：实质上该特性是由Cooley-Tukey算法的本质——分治思想决定的。）

因此在实际使用时，应该严格遵循该性能原则。

另一方面，进一步探究如果数据尺寸不仅仅为质数2、3、5、7的幂次方乘积时，cuFFT的表现如何。

在文档中，可以看到这样一段话：

• 当数据尺寸仅为质数2、3、5、7的幂次方乘积时，cuFFT将调用优化后的 Cooley-Tukey算法计算，此时计算效率最高；
• 当数据尺寸不仅为质数2、3、5、7的幂次方乘积时，还包括其他小于128的质数时，cuFFT将调用普通Cooley-Tukey算法计算；
• 当数据尺寸不仅为质数2到127的幂次方乘积时，cuFFT将调用Bluestein算法计算。需要注意的是，Bluestein算法实现需要比Cooley-Tukey算法实现更多的计算量，且Bluestein算法的精度没有Cooley-Tukey算法高。

在stackoverflow上，有一段关于cuFFT计算效率的描述说的也比较具体：

因此，在使用cuFFT求解傅里叶变换时，需要尽可能按照文档中描述的质数组合，来调整数据的尺寸，以实现算法的最优。

如果涉及到多种尺寸的数据需要同时进行计算dft，注意尽量避免cuFFT调用Bluestein算法计算。对于cuFFT库而言，每当调用Bluestein算法时，都会存在若干隐式的显存开辟和释放的操作，这会严重降低计算速度。

OpenCV中的dft

在OpenCV中，实现了基于CPU的dft运算（cv::dft）和基于CUDA的dft运算（ cv::cuda::dft）。也用到了上述最优数据尺寸的思想。

namespace cv
{

static const int optimalDFTSizeTab[] = {
1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48,
50, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100, 108, 120, 125, 128, 135, 144, 150, 160,
162, 180, 192, 200, 216, 225, 240, 243, 250, 256, 270, 288, 300, 320, 324, 360, 375,
384, 400, 405, 432, 450, 480, 486, 500, 512, 540, 576, 600, 625, 640, 648, 675, 720,
729, 750, 768, 800, 810, 864, 900, 960, 972, 1000, 1024, 1080, 1125, 1152, 1200,
1215, 1250, 1280, 1296, 1350, 1440, 1458, 1500, 1536, 1600, 1620, 1728, 1800, 1875,
1920, 1944, 2000, 2025, 2048, 2160, 2187, 2250, 2304, 2400, 2430, 2500, 2560, 2592,
2700, 2880, 2916, 3000, 3072, 3125, 3200, 3240, 3375, 3456, 3600, 3645, 3750, 3840,
3888, 4000, 4050, 4096, 4320, 4374, 4500, 4608, 4800, 4860, 5000, 5120, 5184, 5400,
5625, 5760, 5832, 6000, 6075, 6144, 6250, 6400, 6480, 6561, 6750, 6912, 7200, 7290,
7500, 7680, 7776, 8000, 8100, 8192, 8640, 8748, 9000, 9216, 9375, 9600, 9720, 10000,
10125, 10240, 10368, 10800, 10935, 11250, 11520, 11664, 12000, 12150, 12288, 12500,
12800, 12960, 13122, 13500, 13824, 14400, 14580, 15000, 15360, 15552, 15625, 16000,
16200, 16384, 16875, 17280, 17496, 18000, 18225, 18432, 18750, 19200, 19440, 19683,
20000, 20250, 20480, 20736, 21600, 21870, 22500, 23040, 23328, 24000, 24300, 24576,
25000, 25600, 25920, 26244, 27000, 27648, 28125, 28800, 29160, 30000, 30375, 30720,
31104, 31250, 32000, 32400, 32768, 32805, 33750, 34560, 34992, 36000, 36450, 36864,
37500, 38400, 38880, 39366, 40000, 40500, 40960, 41472, 43200, 43740, 45000, 46080,
46656, 46875, 48000, 48600, 49152, 50000, 50625, 51200, 51840, 52488, 54000, 54675,
55296, 56250, 57600, 58320, 59049, 60000, 60750, 61440, 62208, 62500, 64000, 64800,
65536, 65610, 67500, 69120, 69984, 72000, 72900, 73728, 75000, 76800, 77760, 78125,
78732, 80000, 81000, 81920, 82944, 84375, 86400, 87480, 90000, 91125, 92160, 93312,
93750, 96000, 97200, 98304, 98415, 100000, 101250, 102400, 103680, 104976, 108000,
···
};
}

int cv::getOptimalDFTSize( int size0 )
{
    int a = 0, b = sizeof(optimalDFTSizeTab)/sizeof(optimalDFTSizeTab[0]) - 1;
    if( (unsigned)size0 >= (unsigned)optimalDFTSizeTab[b] )
        return -1;

    while( a < b )
    {
        int c = (a + b) >> 1;
        if( size0 <= optimalDFTSizeTab[c] )
            b = c;
        else
            a = c+1;
    }

    return optimalDFTSizeTab[b];
}

版权说明

本文为原创文章，独家发布在blog.csdn.net/TracelessLe。未经个人允许不得转载。如需帮助请email至tracelessle@163.com。

参考资料

[1] 快速傅里叶变换 - 维基百科，自由的百科全书
[2] cuFFT :: CUDA Toolkit Documentation
[3] FFTW Home Page
[4] c++ - CUDA FFT - power of two - Stack Overflow
[5] CUDA学习笔记3：CUFFT（CUDA提供了封装好的CUFFT库）的使用例子 - 爱国呐 - 博客园
[6] 快速傅里叶变换 | 晓茵万事通
[7] AutoFFT: A Template-Based FFT Codes Auto-Generation Framework for ARM and X86 CPUs
[8] cuFFT | NVIDIA Developer
[9] FFTW3学习笔记2：FFTW（快速傅里叶变换）中文参考 - 爱国呐 - 博客园
[10] OpenCV: Discrete Fourier Transform
[11] opencv例程解读——dft（离散傅里叶变换）_autocyz-CSDN博客_cv::dft
[12] Understanding profiling (nvprof) output of cuFFT - CUDA / CUDA Programming and Performance - NVIDIA Developer Forums
[13] cuPoisson
[14] 3cHeLoN/cupoisson: CUDA implementation of the 2D fast Poisson solver
[15] cupoisson/poisson.cu at master · 3cHeLoN/cupoisson
[16] Poisson Image Editing
[17] Project 2G: Gradient Domain Editing
[18] opencv/cuda_info.cpp at 68d15fc62edad980f1ffa15ee478438335f39cc3 · opencv/opencv
[19] opencv/dxt.cpp at master · opencv/opencv
[20] OpenCV: Arithm Operations on Matrices
[21] PhDResearchData/cufft.py at 2cf47a41df1e8d9bd9ebbe0da016234bf26cc1d1 · jaisw7/PhDResearchData