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