6

The seminal paper of Cooley and Tukey provides an iterative method for computing the DFT for a sequence of length $N$. Specifically, they mention a method which utilizes the fact that $N$ can be written as $N = r_1 r_2$. The two intermediate results during computations themselves are DFT computations, therefore, their computation in turn can be recursively estimated.

Is there a reference/method of recursive formulation for computing DFT using the $N=r_1 r_2$ factorization ? The versions of FFT I come across use factorizations where $r_1 = \frac{N}{2}$ and $r_2=2$. I would like to know if there are recursive formulations for general factorization of $N$.

curryage
  • 419
  • 4
  • 16
  • @MattL. The versions of FFT I come across use factorizations where $r_1 = \frac{N}{2}$ and $r_2=2$. I would like to know if there are recursive formulations for general factorization of $N$. – curryage Sep 13 '14 at 08:54
  • 1
    The Cooley-Tukey algorithm is not restricted to dividing the transforms into 2 DFTs of length $N/2$. This is just the most popular form. It works for any composite size $N=r_1r_2$. A lot of methods with different factorizations of $N$ are implemented in the FFTW library. You can find more details here. – Matt L. Sep 13 '14 at 09:11
  • 1
    it might be interesting if someone here spelled out how to do it in the case where $r_1 \ne \frac{N}{2}$. – robert bristow-johnson Sep 13 '14 at 16:11
  • 1
    @robertbristow-johnson I finally found the time to do so. – Matt L. Oct 23 '14 at 11:48
  • @MattL., an up vote for that. – robert bristow-johnson Oct 26 '14 at 11:35

2 Answers2

11

The Cooley Tukey method allows general factorizations $N=N_1N_2$. The easiest way to see this is to use index mapping. Consider the length-$N$ DFT of a sequence $x[n]$:

$$X[k]=\sum_{n=0}^{N-1}x[n]W_N^{nk}\tag{1}$$

with $W_N=e^{-j2\pi/N}$. Let's assume that $N$ can be factored as $N=N_1N_2$. Now use the following index mapping

$$n=n_1N_2+n_2,\quad n_1=0,1,\ldots, N_1-1,\quad n_2=0,1,\ldots,N_2-1\\ k=k_1+N_1k_2,\quad k_1=0,1,\ldots, N_1-1,\quad k_2=0,1,\ldots,N_2-1\tag{2}$$

Plugging (2) into (1) gives

$$X[k_1+N_1k_2]=\sum_{n_2=0}^{N_2-1}\sum_{n_1=0}^{N_1-1}x[n_1N_2+n_2]W_N^{(n_1N_2+n_2)(k_1+k_2N_1)}\tag{3}$$

The complex exponentials in (3) can be simplified as follows:

$$W_N^{(n_1N_2+n_2)(k_1+k_2N_1)}=W_N^{N_2n_1k_1}W_N^{N_1N_2n_1k_2}W_N^{n_2k_1}W_N^{N_1n_2k_2}=\\= W_{N_1}^{n_1k_1}\cdot 1\cdot W_{N}^{n_2k_1}W_{N_2}^{n_2k_2}\tag{4}$$

because $W_N^{N_1}=W_{N_2}$, $W_N^{N_2}=W_{N_1}$, and $W_N^{N_1N_2}=W_N^N=1$. Plugging (4) back into (3) and rearranging terms gives

$$X[k_1+N_1k_2]=\sum_{n_2=0}^{N_2-1}\left[W_N^{n_2k_1}\left(\sum_{n_1=0}^{N_1-1}x[n_1N_2+n_2]W_{N_1}^{n_1k_1}\right)\right]W_{N_2}^{n_2k_2}\tag{5}$$

Equation (5) can be interpreted as follows: the result can be computed by computing $N_2$ length-$N_1$ DFTs of subsampled versions of the input sequence $x[n]$ (this is the inner sum in (5)), then multiplying the result by twiddle factors $W_N^{n_2k_1}$, and finally computing $N_1$ length-$N_2$ DFTs (the outer sum in (5)). This procedure is most easily visualized using matrices:

  1. read the data row-wise into a $N_1\times N_2$ matrix
  2. compute DFTs over all $N_2$ columns (this can be done in-place)
  3. multiply the result element-wise with twiddle factors
  4. compute DFTs over all $N_1$ rows (in-place)
  5. read out the data column-wise to get the final length-$N$ result

This little octave/Matlab script illustrates the procedure:

N = 1536;          % non-power of 2 FFT length used in LTE
N1 = 512; N2 = 3;  % 512 x 3 = 1536
x = randn(N,1);    % some input signal

% generate N1 x N2 matrix of input data, reading them row wise
X = reshape(x,[N2,N1]).';

% DFT over each column
X = fft(X);

% Multiplication with twiddle factors
X = X.*exp(-1i*2*pi/N*(0:N1-1)'*(0:N2-1));

% DFT over each row
X = (fft(X.')).';

% read out column wise
X = X(:);
Matt L.
  • 89,963
  • 9
  • 79
  • 179
  • 1
    Excellent answer. I just want to also add a link to the derivation I found helpful in making the mixed radix kissfft http://sourceforge.net/projects/kissfft/files/derivation/ – Mark Borgerding Oct 23 '14 at 13:40
2

Matt's answer is excellent. For anyone wanting to do the same in Python, I thought I'd add the code (with a slight modification to use a sine wave input, which I prefer for debugging):

import numpy as np

N=1536
N1=512
N2=3                    # 512 x 3 = 1536
nCycles = 2             # Number of cycles in source frame


time=np.arange(0,N )    # some input signal
x=np.sin(nCycles * time * 2. * np.pi / N)


generate N1 x N2 matrix of input data, reading them row wise


x=np.reshape(x,(N1,N2))


DFT over each column


X=np.fft.fft(x,axis=0)


Element-wise multiplication with twiddle factors


X=np.multiply(X,np.exp(- 1j * 2. * np.pi / N * (np.arange(0,N1).T.reshape(N1,1) @ np.arange(0,N2).reshape(1,N2))))


DFT over each row


X=np.fft.fft(X,axis=1)


read out column wise


X=np.ravel(X.T)
print(np.round(X,4))
Johned
  • 31
  • 1