Finding Suitable FFT Numbers

sfft5

Many high quality FFT libraries support fast implementations for small primes. (2,3,5,7). I’ll use CUDA fft here as an example.

When you have say  an array of 18000 samples d_data, you might not want to do this:

nfft = 32768; 
cudafftPlan1d(&m_plan, nfft, CUFFT_C2C, 1); 
cufftExecC2C(m_plan, d_data, d_data, CUFFT_FORWARD);

Especially if you are doing many, many iterations of this.

Instead, you can use a method to find the closest fft number to 18000, that is still made of of small primes (2,3,5,7).

With this findNFFT method, nfft = 18432, which is

2^11 * 3 * 3

In this way, you process less samples, which saves on time. Any gain in speed by using a true radix-2 (32768) number is offset by the number of samples being processed, which is now much less. In one of my projects, this helped to cut the processing time by half.

Here’s the code for the findNFFT method – not the prettiest code but it does the trick.

/**
 * @brief Returns a vector containing the prime factors of n
 *
 * @param [in] The number to find the prime factors for
 * @return
 */
std::vector<int> primeFactors(int n) {
    std::vector<int> vec;

    while (n % 2 == 0) {
        vec.push_back(2);
        n /= 2;
    }

    for (int i = 3; i <= sqrt(n); i += 2) {
        while (n % i == 0) {
            vec.push_back(i);
            n /= i;
        }
    }

    if (n > 2)
        vec.push_back(n);

// std::cout << "Prime factors:" << std::endl;
// for (int j=0; j < vec.size(); j++)
// {
// printf("%d ", vec[j]);
// }
// printf("\n");
    return vec;
}

/**
 * @brief Used to find the appropriate fft integer for the input n
 * This uses the "formula" (N + D - 1)/D * D
 * Criteria: Output nfft should be a factor of 2,3,5
 *
 * @param [in] Integer to find nfft for
 */
int findNFFT(int n) {
    std::vector<int> ansPrimes;
    std::vector<int> firstPrimes;

    int d = 0;

    do {
        if (n > 2048) d = 512;
        else if (n > 1024) d = 256;
        else if (n > 128) d = 64;
        else if (n > 32) d = 32;
        else if (n > 8) d = 8;
        else d = 2;

        int fn = (n + d - 1) / d * d;
        firstPrimes = primeFactors(fn);

        for (int i = 0; i < firstPrimes.size(); i++) {
            if (firstPrimes[i] == 2 || firstPrimes[i] == 3 || firstPrimes[i] == 5) {
                ansPrimes.push_back(firstPrimes[i]);
                firstPrimes.erase(firstPrimes.begin() + i);
                i -= 1;
            }
        }

        int newN = 1;
        if (firstPrimes.size() > 0) {
            for (int i = 0; i < firstPrimes.size(); i++)
                newN *= firstPrimes[i];
        }

        n = newN;
        firstPrimes = {};

    } while (n != 1); // if n == 1 means that firstPrimes

    int ans = 1;
    for (int i = 0; i < ansPrimes.size(); i++)
        ans *= ansPrimes[i];

    return ans;
}

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s