FFT with powers of 3

Question

Classic Fast Fourier Transfrom (FFT) works fine, when $n$ is power of 2. How to generalize FFT procedure when $n$ is power of 3? Is it possible to easily modify the algorithm and preserve its correctness?

Thanks in advance.

Answer 1

FFT is defined for every decomposition to prime factors, i.e., if $$ n=p_1^{r_1}\cdots p_k^{r_k}, $$ then the FFT of an $n$-vector is definable in $r_1+\cdots+r_k$ steps:

Step 1. We create a $p_1$-vector,

Step 2. We create a $p_1^2$-vector,...,

Step $r_1$. We create a $p_1^{r_1}$-vector,...,

Step $r_1\!+\!1$. We create a $p_1^{r_1}p_2$-vector, etc.

Answer 2

Let me try to explain my view point for the recursive step for the FFT computation. First define a few things.

Let $F_N=[\omega_n^{2\pi ikl/N}]_{k,l\leq N}$, be the FFT matrix of order $N$ and $\omega_N = e^{2\pi i /N}$ is an $N$'th root of unity.

Let $A\in{\mathbb R}^{m_1\times n_1}$ and $B\in{\mathbb R}^{m_1\times n_1}$ then the Kronecker product $A\otimes B$ and box product $A\boxtimes B$ are defined by $$ (A\otimes B)_{(i-1)m_2+j,(k-1)n_2+l}=a_{ik}b_{jl}=(A\otimes B)_{(i,j),(k,l)} $$ and $$ (A\boxtimes B)_{(i-1)m_2+j,(k-1)n_1+l}=a_{il}b_{jk}=(A\otimes B)_{(i,j),(k,l)} $$ These two operations are useful because $(A\otimes B)\mathrm{vec}(X)=\mathrm{vec}(AXB)$ and $(A\boxtimes B)\mathrm{vec}(X)=\mathrm{vec}(AX^\top B)$.

Now define $$ V_{mk}(\alpha)=\begin{pmatrix} 1 & 1 & 1 & \cdots & 1\\ 1 & \alpha & \alpha^2 & \cdots & \alpha^{k-1}\\ 1 & \alpha^2 & \alpha^4 & \cdots & \alpha^{2(k-1)}\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \alpha^{m-1} & \alpha^{2(m-1)} & \cdots & \alpha^{(m-1)(k-1)} \end{pmatrix} $$ then if $N=km$ $$ F_N=(F_k\otimes I_m)\mathrm{diag}(V_{mk}(\omega_N))(I_k\boxtimes F_m) $$ To compute the FFT of a signal $x$ we simply need to compute $F_N x$. Now reshape $x$ into a matrix $X$ such that $x=\mathrm{vec}(X)$ then we can use the rules for the Kronecker and box product to compute the matrix vector product efficiently: $$ F_N \mathrm{vec}(X) = \mathrm{vec}((V_{mk}(\omega_N))\circ (F_m X^\top))F_k^\top), $$ where $\circ$ denotes element-wise multiplication. Instead of taking $O(N^2)=O(k^2m^2)$ operations this expression does the job in $O(km(k+m))$ operations. Now simply repeat the procedure to factorize $F_k$ and $F_m$ recursively.

Answer 3

I don't know what you want to do with FFT, so this might not be relevant for you,

A nice generalization of FFT to arbitrary size, not just powers of two, is the truncated Fourier transform, introduced by Joris van der Hoeven which behaves well for all size, and smoothes the jumps of the FFT. It applies very well to multiplication algorithms (polynomial, integers) based on Fourier transform.

Link to the paper : http://www.texmacs.org/joris/issac04/issac04.pdf

Answer 4

Disappointed that it took so long to properly answer this. The question is simple and the algorithm it is asking to implement isn't that hard. I just implemented it in Python. Pasting the code below. Will add more explanation later, but reading chapter 30 of the book "Introduction to algorithms" by CLRS should make things clear.

def fft_3(x):
    """
    Perform a fast Fourier transform
    on the given array of complex numbers.
    """ 
    n = len(x)
    if n == 1: 
        return x
    # The complex roots of unity.
    omega_n = np.exp(1j*2*np.pi/n)
    omega = 1+0j
    # Extact the even indices from array.
    y_0 = fft_3(x[::3]) 
    # Extract the odd indices from array.
    y_1 = fft_3(x[1::3])
    y_2 = fft_3(x[2::3]) 
    # The result array.
    y = [0]*n
    # The step that merges the two inputs into the
    # final answer. This step takes O(n) time since
    # it is just a single loop.
    for k in range(n//3):
        y[k] = y_0[k]+omega*y_1[k]+omega**2*y_2[k]
        y[k+n//3] = y_0[k]+np.exp(2*np.pi*1j/3)*omega*y_1[k]+np.exp(4*np.pi*1j/3)*omega**2*y_2[k]
        y[k+(2*n)//3] = y_0[k]+np.exp(4*np.pi*1j/3)*omega*y_1[k]+np.exp(2*np.pi*1j/3)*omega**2*y_2[k]
        omega = omega*omega_n
return y

Just like the original FFT in most books only works on powers of $2$ and will produce incorrect results otherwise, this version will only work on powers of $3$. A basic sanity check:

import numpy as np
fft_3([1,2,3,4,5,6,7,8,9])
np.fft.fft([1,2,3,4,5,6,7,8,9])

The sign of the complex part of my algorithm is the opposite of the numpy version, but that is probably an artifact of convention and can be switched trivially.