6
$\begingroup$

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$.

$\endgroup$
5
  • $\begingroup$ @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$. $\endgroup$
    – curryage
    Sep 13, 2014 at 8:54
  • 1
    $\begingroup$ 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. $\endgroup$
    – Matt L.
    Sep 13, 2014 at 9:11
  • 1
    $\begingroup$ it might be interesting if someone here spelled out how to do it in the case where $r_1 \ne \frac{N}{2}$. $\endgroup$ Sep 13, 2014 at 16:11
  • 1
    $\begingroup$ @robertbristow-johnson I finally found the time to do so. $\endgroup$
    – Matt L.
    Oct 23, 2014 at 11:48
  • $\begingroup$ @MattL., an up vote for that. $\endgroup$ Oct 26, 2014 at 11:35

2 Answers 2

11
$\begingroup$

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(:);
$\endgroup$
1
2
$\begingroup$

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))
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.