FFTC1: One-Dimensional Complex Fourier Transforms

Whether you prove them to yourself or take my word for them the formulas given above should be enough for you to follow what the routine fftc1 does. The first section takes an array of $N$ points and rearranges them into bit-reversed order. This is done by noting that if you know the bit-reversal $BR_b$ of a binary number $b$ it is easy to find the bit-reversal $BR_{b+1}$ of the number $b+1$. To add $1$ to $b$ you start at the right and replace every $1$ with a $0$ until you reach the first $0$, which you replace with a $1$. So to go from $BR_b$ to $BR_{b+1}$ you do the bit-reverse of this process, starting at the left and replacing every $1$ with a $0$ until you reach the first $0$, which you replace with a $1$.

The next part of the routine implements the DL formula on successively larger blocks, beginning with combining $N=1$ transforms into $N=2$ transforms. Each time the routine calculates a transform of size $trans\_size$ the first $trans\_size/2$ elements represent $F^e_b$ and the next half represent $F^o_b$. These are then used to calculate $F_b$ and $F_{b+trans\_size/2}$. Note that in the DL formula $F^e_{b+trans\_size/2}=F^e_b$, $F^o_{b+trans\_size/2}=F^o_b$, and $W^{b+trans\_size/2} = - W^b$, so the only difference between the formulas for $F_b$ and $F_{b+trans\_size/2}$ is the sign of $W^b$. Finally, the inner loop over the variable $trans$ iterates over all the transforms of size $trans\_size$ in the entire array, so if you are transforming an $N=128$ array of data you first have to calculate $64$ $trans\_size=2$ transforms, then $32$ $trans\_size=4$ transforms, and so on.

Note that the formulas for the real and imaginary parts of $wb$ are nothing but $\cos(2 \pi b/N)$ and $\sin(2 \pi b/N)$. They could have simply been calculated directly but for the sake of efficiency they are calculated instead using a recursion formula each time $b$ is incremented. This formula can be verified simply using the addition rules for sines and cosines.

Note that the end result of this process is a transform with the results in sequential order from $b=0$ to $b=N-1$, which by the periodicity of $F_b$ is exactly the storage scheme described in section 2.3.2.

The reverse transform is exactly identical except the sign of the exponential is reversed and the final output is divided by $N$.

The argument $skip$ is used in case the elements of the input array are not sequential but evenly spaced in some larger array. To do a one-dimensional transform you would almost always want to set $skip$ to one, but this option allows fftc1 to be used by fftcn, which does transforms of multi-dimensional arrays. For all but the last array index the one-dimensional transforms used by fftcn are not of sequential data.