« FFT Zero-Padding...
FFT Zero-Padding... »

FFT Zero-Padding Strategies - Small Prime Factors

Steve Eddins
Posted on Thursday, May 2, 2024

This post is part of the FFT Zero-Padding Strategies series.

Power-of-two padding

In my previous post, I wrote about using FFTs to compute a full convolution. If the vector x has $K$ points and the vector h has $L$ points, then the convolution of x and h has $K+L-1$ points. FFTs can be used to compute this convolution, but only if x and h are zero-padded to have at least $K+L-1$ points in the FFT computation.

Early in my career, the commonly available FFT algorithms, which were usually written in Fortran, were limited to computing FFT with transforms that were powers of two. Because of this limitation, using FFTs to compute convolution required zero-padding to the smallest power of two that was greater than or equal to $K+L-1$ . Here’s how that might look in MATLAB code.

x = [1 -1 0 4];
h = [1 0 -1];
K = length(x);
L = length(h);

N = K + L - 1
N = 6
Np = 2^nextpow2(N)
Np = 8
X = fft(x,Np);   % Computes Np-point zero-padded FFT
H = fft(h,Np);   % Computes Np-point zero-padded FFT

Y = ifft(X .* H);
y = Y(1:N)
y = 1x6    
1.0000   -1.0000   -1.0000    5.0000         0   -4.0000

Small-primes padding

More than 20 years ago now, I integrated the FFTW library into MATLAB. FFTW, or “Fastest Fourier Transform in the West,” won the 1999 J. H. Wilkinson Prize for Numerical Software. I spent a long time deep in the FFTW documentation at the time, and this statement caught my attention: “The standard FFTW distribution works most efficiently for arrays whose size can be factored into small primes (2, 3, 5, and 7), and otherwise it uses a slower general-purpose routine.”

This observation suggests the interesting possibility of using transform lengths are powers of 2, 3, 5, and 7, instead of using only transform lengths that are powers of 2.

Here is a function that computes a transform length based on this idea.

function np = transformLength(n)
    np = n;

    while true
        % Divide n evenly, as many times as possible, by the factors 2, 3,
        % 5, and 7.
        r = np;
        for p = [2 3 5 7]
            while (r > 1) && (mod(r, p) == 0)
                r = r / p;

        if r == 1
            % If the result after the above divisions is 1, then we have found
            % the desired number, so break out of the loop.
            % np has one or more prime factors greater than 7, so try the
            % next integer.
            np = np + 1;

Suppose n is 37, a prime number. The next power of 2 is 64. The next number whose largest prime factor 7 or less, on the other hand, is 40.

n = 37;
np = transformLength(n)
np = 40
ans = 1x4    
     2     2     2     5

2-D performance comparisons

To see how different padding strategies affect performance, I will set up some different FFT-based 2-D convolution functions and measure their execution time for different input sizes. Here is function that zero-pads only to the minimum amount required, $K+L-1$ .

function C = conv2_fft_minpad(A,B)
    [K1,K2] = size(A);
    [L1,L2] = size(B);
    N1 = K1 + L1 - 1;
    N2 = K2 + L2 - 1;

    C = ifft2(fft2(A,N1,N2) .* fft2(B,N1,N2));

Let’s measure the time required to convolve a 1170x1170 input with a 101x101 input. The output will be 1271x1271, and 1271 has prime factors 31 and 41.

A = rand(1170,1170);
B = rand(101,101);
f_minpad = @() conv2_fft_minpad(A,B);
ans = 0.0329

Here is a function that zero-pads to the next power of 2. (The next power of 2 above 1271 is 2048.)

function C = conv2_fft_pow2pad(A,B)
    [K1,K2] = size(A);
    [L1,L2] = size(B);
    N1 = K1 + L1 - 1;
    N2 = K2 + L2 - 1;
    N1p = 2^nextpow2(N1);
    N2p = 2^nextpow2(N2);

    C = ifft2(fft2(A,N1p,N2p) .* fft2(B,N1p,N2p));
    C = C(1:N1,1:N2);

Repeat the timing experiment using this second padding method:

f_pow2pad = @() conv2_fft_pow2pad(A,B);
ans = 0.0553

So, power-of-2 padding is not helping in this case. The fact that power-of-2 padding is slower than using a transform length with relatively high prime factors represents a sea change from the way FFT computations in MATLAB performed prior to 2004, which is when MATLAB started using FFTW.

Can small-primes padding do better?

function C = conv2_fft_smallprimespad(A,B)
    [K1,K2] = size(A);
    [L1,L2] = size(B);
    N1 = K1 + L1 - 1;
    N2 = K2 + L2 - 1;
    N1p = transformLength(N1);
    N2p = transformLength(N2);

    C = ifft2(fft2(A,N1p,N2p) .* fft2(B,N1p,N2p));
    C = C(1:N1,1:N2);

Run the timing again with the small-primes pad method:

f_smallprimespad = @() conv2_fft_smallprimespad(A,B);
ans = 0.0230

Yes, the convolution function using the small-primes method runs about 45% faster for this case.

Now let’s compare the results for a range of sizes.

nn = 1100:1200;
t_minpad = zeros(size(nn));
t_pow2pad = zeros(size(nn));
t_smallprimespad = zeros(size(nn));
for k = 1:length(nn)
    n = nn(k);
    A = rand(n,n);
    t_minpad(k) = timeit(@() conv2_fft_minpad(A,B));

    t_pow2pad(k) = timeit(@() conv2_fft_pow2pad(A,B));

    t_smallprimespad(k) = timeit(@() conv2_fft_smallprimespad(A,B));

Plot the results.

ax = gca;
ax.YLim(1) = 0;
legend(["min padding","power-of-2","small-primes"],...
title("Execution time, nxn convolution with 101x101")
ylabel("time (s)")
grid on


While this is admittedly not anything like a comprehensive performance study, it does suggest that:

In my next post in this series, I plan to investigate a possible way to speed up the computation of the small-primes transform length.

This series

  1. FFT Zero-Padding Strategies - Introduction
  2. FFT Zero-Padding Strategies - Small Prime Factors
  3. FFT Zero-Padding Strategies - Computing Transform Length