FFT butterfly input index

I am trying to determine a "simple" way to compute which inputs of a FFT need to "butterfly" together for its various stages. I am looking at a diagram like this:

In that diagram of a $16$ point FFT,

• For stage 1: 0 butterflies with 8, 4 butterflies with 12, and so on.
• For stage 2: 0 butterflies with 4, 8 butterflies with 12, and so on. Of course stage 2's inputs are stage 1's outputs.

My hope is that for each stage of the FFT, I can just have a simple counter going from $0$ to $N-1$ (FFT length $N$), and I can do very inexpensive bit manipulation to form the correct indices in the correct order.

Stage 1's pattern is fairly obvious: Have a 4 bit counter going from 0 to 15, and completely reverse the bits in the counter:

• $0000 (0)$ with $0001\longleftrightarrow 1000(8)$
• $0010\longleftrightarrow 0100(4)$ with $0011\longleftrightarrow 1100 (12)$
• and so on.

I believe stage 2 is like this: bit reverse the counter, then swap the 2 leftmost bits.

Is there a generic formula for any stage of any FFT?

• What are the inputs and desired outputs in your desired formula? For example, inputs: $N=16$, $s=1$, and the output should be $[0, 8, 4, 12, 2, 10, 6, 14, 1, 9, 5, 13, 3, 11, 7, 15]$? or for input $N=16$, $s=2$, the output should be $[0, 4, 2, 6, 1, 5, 3, 7, 8, 12, 10, 14, 9, 13, 11, 15]$? (where the pairs are given sequentially) Is that what you are looking for? – msm Dec 2 '16 at 6:36
• yes! that is exactly what I am looking for. – user2913869 Dec 2 '16 at 14:39

Assume $N=16$. so the counter $c$ chooses values between $0$ to $15$ and there are $s=1,2,\cdots,\log2(N)$ number of stages. Based on your comment, the input-output relationship for stage $s$ can be given by:

• Represent the counter in binary form in $\log_2 N$ bits (here 4 bits).

• Circularly shift the counter $s-1$ times to the right.

• Bit reverse the result.

• Represent in decimal.

The following Matlab code can give the output in vector form as explained in the comment:

N = 16;
stage = 1              %Do it here for stage #1
c = 0:N-1;
c_bin = de2bi(c);
butterfly_idx = bi2de(fliplr(circshift(c_bin',stage-1)'))


Example for N=16 and stage=2:

c_bin =

0     0     0     0
1     0     0     0
0     1     0     0
1     1     0     0
0     0     1     0
1     0     1     0
0     1     1     0
1     1     1     0
0     0     0     1
1     0     0     1
0     1     0     1
1     1     0     1
0     0     1     1
1     0     1     1
0     1     1     1
1     1     1     1

fliplr(circshift(c_bin',1)')

ans =

0     0     0     0
0     0     1     0
0     1     0     0
0     1     1     0
1     0     0     0
1     0     1     0
1     1     0     0
1     1     1     0
0     0     0     1
0     0     1     1
0     1     0     1
0     1     1     1
1     0     0     1
1     0     1     1
1     1     0     1
1     1     1     1

bi2de(fliplr(circshift(c_bin',1)'))

ans =

0
4
2
6
1
5
3
7
8
12
10
14
9
13
11
15


which is the expected result for the second stage.

• I tested this, and it seems almost correct. The one thing that threw me off is that my de2bi was returning binary numbers in bit-reversed order (LSb was on the left), so when you did a right shift, that would normally equate to a left shift on the non reversed binary number. So I think "times to the right" should say "times to the left" – user2913869 Dec 7 '16 at 18:40
• It is actually a little tricky. First, you can control the position of the LSB by adding right-msb or left-msb a second argument to de2bi. Second, the matrix is transposed before it is given as an argument and the result is again transposed. I have tested this on Matlab 2010 and it works as expected. But I am not sure how this works in new versions. In the setting I explains wich was assumed to be default, LSB should be on the right (as you confirmed) and after each circular shift the LSB should go to the right and the MSB should be fed to the left-most position. I will add an example... – msm Dec 7 '16 at 21:39

So this FFT code from loooong ago is about as simple as it can be that implements both the Decimation-in-Frequency (the FFT) and Decimation-in-Time (the iFFT). i think you'll be able to find a piece of code that matches your N=16 FFT diagram.

/*
A set of utility programs to compute the Fast Fourier Transform (FFT):

N-1
X[k] =     SUM { x[n]exp(-j2 pi nk/N) }
n=0

and inverse Fast Fourier Transform (iFFT):

N-1
x[n] = 1/N SUM { X[k]exp(+j2 pi nk/N) }
k=0

To speed things up, multiplication by 1 and j are avoided.  The input, x[], is an array of
complex numbers (pairs of doubles) of length N = 2^n.  The calling program supplies
n = log2(N) not the array length, N.  The output is placed in BIT REVERSED order in x[].
Call bit_reverse(x, n) to swap back to normal order, if needed. However, iFFT(X, n, stbl)
requires its input, X[], to be in bit reversed order making bit reversing unnecessary in
some cases, such as convolution.  stbl is an array of doubles of length N/4 containing the
sine function from 0 to pi/2 used to compute the FFT.  Call sin_table(stbl, n) ONLY ONCE
before either FFT(x, n, stbl) or iFFT(X, n, stbl) to set up a sin table for FFT computation.

Written ca. 1985 in THINK C by Robert Bristow-Johnson.
*/

#define HALFPI 1.570796326794897
#define PI     3.141592653589793
#define TWOPI  6.283185307179586

// #include "complex.h"

typedef struct {
double real;
double imag;
} complex;

#define Re(z) (z).real
#define Im(z) (z).imag

void FFT(complex *x, int n, double *stbl)
{
long size;
register long length, step, stepsize, end;
register complex *top, *bottom;                                       /* top & bottom of FFT butterfly */
complex temp;

size = 1L<<(n-2);
end = (long)x + 4*sizeof(temp)*size;

length = size;
stepsize = 1;
while ( length >= 1)
{
top = x;
while ((long)top < end)
{
bottom = top + 2*length;

Re(temp) = Re(*top) - Re(*bottom);                    /* butterfly: twiddle = 1 */
Im(temp) = Im(*top) - Im(*bottom);
Re(*top) += Re(*bottom);
Im(*top) += Im(*bottom);
*bottom = temp;
top++;
bottom++;

for (step = stepsize; step < size; step += stepsize)
{
Re(temp) = Re(*top) - Re(*bottom);            /* butterfly: twiddle = exp(-j theta) */
Im(temp) = Im(*top) - Im(*bottom);
Re(*top) += Re(*bottom);
Im(*top) += Im(*bottom);
Re(*bottom) = Re(temp)*stbl[size - step] + Im(temp)*stbl[step];
Im(*bottom) = Im(temp)*stbl[size - step] - Re(temp)*stbl[step];
top++;
bottom++;
}

Re(temp) = Im(*top) - Im(*bottom);                    /* butterfly: twiddle = -j */
Im(temp) = Re(*bottom) - Re(*top);
Re(*top) += Re(*bottom);
Im(*top) += Im(*bottom);
*bottom = temp;
top++;
bottom++;

for (step = stepsize; step < size; step += stepsize)
{
Re(temp) = Im(*top) - Im(*bottom);            /* butterfly: twiddle = -j*exp(-j theta) */
Im(temp) = Re(*bottom) - Re(*top);
Re(*top) += Re(*bottom);
Im(*top) += Im(*bottom);
Re(*bottom) = Re(temp)*stbl[size - step] + Im(temp)*stbl[step];
Im(*bottom) = Im(temp)*stbl[size - step] - Re(temp)*stbl[step];
top++;
bottom++;
}
top = bottom;
}
length >>= 1;
stepsize <<= 1;
}

top = x;
bottom = x + 1;
while ((long)top <  end)
{
Re(temp) = Re(*top) - Re(*bottom);                            /* butterfly: twiddle = 1 */
Im(temp) = Im(*top) - Im(*bottom);
Re(*top) += Re(*bottom);
Im(*top) += Im(*bottom);
*bottom = temp;
top += 2;
bottom += 2;
}
}

void iFFT(complex *X, int n, double *stbl)
{
long size;
register long length, step, stepsize, end;
double scale;
register complex *top, *bottom;                                       /* top & bottom of FFT butterfly */
complex temp;

size = 1L<<(n-2);
end = (long)X + 4*sizeof(temp)*size;

scale = 0.25/size;
top = X;
bottom = X + 1;
while ((long)top <  end)
{
Re(temp) = (Re(*top) - Re(*bottom))*scale;                    /* butterfly: twiddle = 1/N */
Im(temp) = (Im(*top) - Im(*bottom))*scale;
Re(*top) = (Re(*top) + Re(*bottom))*scale;
Im(*top) = (Im(*top) + Im(*bottom))*scale;
*bottom = temp;
top += 2;
bottom += 2;
}

length = 1;
stepsize = size;
while ( stepsize >= 1)
{
top = X;
while ((long)top < end)
{
bottom = top + 2*length;

temp = *bottom;                                       /* butterfly: twiddle = 1 */
Re(*bottom) = Re(*top) - Re(temp);
Im(*bottom) = Im(*top) - Im(temp);
Re(*top) += Re(temp);
Im(*top) += Im(temp);
top++;
bottom++;

for (step = stepsize; step < size; step += stepsize)
{                                             /* butterfly: twiddle = exp(+j theta) */
Re(temp) = Re(*bottom)*stbl[size - step] - Im(*bottom)*stbl[step];
Im(temp) = Im(*bottom)*stbl[size - step] + Re(*bottom)*stbl[step];
Re(*bottom) = Re(*top) - Re(temp);
Im(*bottom) = Im(*top) - Im(temp);
Re(*top) += Re(temp);
Im(*top) += Im(temp);
top++;
bottom++;
}

Re(temp) = -Im(*bottom);                              /* butterfly: twiddle = +j */
Im(temp) = Re(*bottom);
Re(*bottom) = Re(*top) - Re(temp);
Im(*bottom) = Im(*top) - Im(temp);
Re(*top) += Re(temp);
Im(*top) += Im(temp);
top++;
bottom++;

for (step = stepsize; step < size; step += stepsize)
{                                             /* butterfly: twiddle = +j*exp(+j theta) */
Re(temp) = -Im(*bottom)*stbl[size - step] - Re(*bottom)*stbl[step];
Im(temp) = Re(*bottom)*stbl[size - step] - Im(*bottom)*stbl[step];
Re(*bottom) = Re(*top) - Re(temp);
Im(*bottom) = Im(*top) - Im(temp);
Re(*top) += Re(temp);
Im(*top) += Im(temp);
top++;
bottom++;
}
top = bottom;
}
length <<= 1;
stepsize >>= 1;
}
}

void sin_table(double *stbl, int n)
{
register long size, i;
double theta;

size = 1L<<(n-2);
theta = HALFPI/size;

for (i = 0; i < size; i++)
{
stbl[i] = sin(theta*i);
}
}

void bit_reverse(register complex *x, int n)
{
complex temp;
register long k, i, r, size, count;

size = (1L<<n) - 1L;
for (k = 1L; k < size; k++)
{
i = k;
r = 0;
for (count = n; count > 0; count--)
{
r <<= 1;
r += (i & 0x00000001L);
i >>= 1;
}
if (r > k)
{
temp = x[r];
x[r] = x[k];
x[k] = temp;
}
}
}

• This answer is unclear to me. Can you show specifically in the code where it addresses my original question? – user2913869 Dec 2 '16 at 16:14
• okay, so my iFFT() code is Decimation-in-Time, what your FFT diagram is. "Stage 1" is qualitatively different than the other stages in that there is not a butterfly twiddle that is off by a factor of $j = e^{j \pi/4}$. your "Stage 1" is where my comment /* butterfly: twiddle = 1/N */ is. the dark lines in your later stages are where my comment /* butterfly: twiddle = exp(+j theta) */ is and the butterflies where your grey lines are are where my /* butterfly: twiddle = +j*exp(+j theta) */ is. – robert bristow-johnson Dec 2 '16 at 18:07
• but in both cases where /* butterfly: twiddle = exp(+j theta) */ and /* butterfly: twiddle = +j*exp(+j theta) */ when theta is zero, my old code saves on multiplication with twiddle = 1 and twiddle = j`. this was because, back 30 years ago, multiplication cost much more than addition. so maybe i should take that code out. – robert bristow-johnson Dec 2 '16 at 18:11