This article describes the design, implementation and testing of a C# program implementing a fundamental operation in the computation of the Fast Fourier Transform.
Introduction
Given a set of N complex-valued data points (usually samples of a signal at discrete times), denoted by hk for
, their discrete Fourier transform (DFT) is defined as:

(A detailed description of the DFT is beyond the scope of this article. The reader is referred to the book Numerical Recipes in C++, The Art of Scientific Computing by W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery, Cambridge University Press, Second Edition, 2002 (hereafter abbreviated as “The Reference”), where the equation above is numbered 12.1.7).
The DFT can be readily implemented on any modern programming language, but its execution can be shown to be quadratic in N, denoted by
. There are several algorithms to improve the execution time to compute the DFT.
The fast Fourier transform (FFT) computes the DFT in
time, and is defined as follows. Defining the complex number W as:

the equation for the DFT can be written as:

(These are, respectively, equations 12.2.1 and 12.2.2 in The Reference.)
One of the implementations of the FFT, as described by the previous equation, was “re-discovered” by G.C. Danielson and C. Lanczos in 1942 (“Some improvements in practical Fourier analysis and their application to X-ray scattering from liquids,” J. Franklin Inst. 233, 365-380 and 435-452). They showed that a DFT of power-of-two length N can be re-written as the sum of two DFTs, each of length N/2. One DFT is computed from the even-numbered data points; the other is computed from the odd-numbered points. Their derivation of the algorithm is as follows (see The Reference, p. 510, 12.2.3).

In the derivation, W is the same complex constant defined before,
is the kth element of the DFT of length N/2 for the even-indexed original data points
, and
is the kth element of the DFT of length N/2 for the odd-indexed data points.
Implementation of the FFT
The Danielson-Lanczos algorithm can be implemented recursively. The base case of the recursion occurs when there is just one data point, whose transform is itself. For every pattern of
e’s and o’s, there is a one point transform that is just one of the input data points
for some n:

The preceding one-point transform is independent of k because it is periodic in k with period 1. To determine which value of n corresponds to a given pattern of e’s and o’s in the last equation, reverse the pattern of e’s and o’s, let e = 0 and o = 1, and the result is the value of n in binary. This process is called bit reversal and it is the preparatory step before applying the Daniellson-Lanczos algorithm to the data points. The Reference (p. 511) shows a graphical representation of the process:

The figure depicts the re-ordering of an array of complex numbers by bit reversal, either using two arrays (left) or in place (right). Anticipating the C++ implementation of the FFT in The Reference, N complex numbers are stored in an array of size 2N. For each number, the real part is stored in one element and the imaginary part in the next element. Thus the first complex number,
, would be stored in array elements 0 and 1.
The bit-reversal process is illustrated with arrows on the left side of the figure above. For example, the binary index 001 is mapped to index 100, the index 011 is mapped to index 110, and so on. Symmetric indices are mapped to themselves. As will be seen later, the mappings are intended to indicate swap operations. That is, if integer-valued index j is mapped to index k, the array elements at those indices must be interchanged, as depicted on the right side of the figure above. However, there is a problem. Observe that in the two cases shown, the real part of a complex number is interchanged with the imaginary part of another number. Clearly, the swaps are not being done properly.
The Reference goes on (pp. 511-512) explaining what the FFT does and how the complex data points are stored in a real-valued C++ array. If the number of complex data points is nn, the actual length of the real-valued array is 2 times nn (data[ 0 .. 2*nn – 1), with each complex number occupying two consecutive locations: data[ 0 ] is the real part of
, data[ 1 ] is the imaginary part of
, and so on.
The FFT routine overwites the data array with the DFT complex values. The complete C++ code implementing the FFT is available in the public domain, and can be obtained from (among other sites):
Of the code downloaded by the author, function four1
is the one that implements the FFT:
#include <cmath>
#include "nr.h"
using namespace std;
void NR::four1( Vec_IO_DP &data, const int isign )
{
int n, mmax, m, j, istep, i;
DP wtemp, wr, wpr, wpi, wi, theta, tempr, tempi;
int nn = data.size() / 2;
n = nn << 1;
j = 1;
for ( i = 1; i < n; i += 2 ) {
if ( j > i )
{
SWAP( data[ j - 1 ], data[ i - 1 ] ); SWAP( data[ j ], data[ i ] );
}
m = nn;
while ( m >= 2 && j > m )
{
j -= m;
m >>= 1;
}
j += m;
}
mmax = 2;
while ( n > mmax ) { istep = mmax << 1;
theta = isign * ( 6.28318530717959 / mmax ); wtemp = sin( 0.5 * theta ); wpr = -2.0 * wtemp * wtemp;
wpi = sin( theta );
wr = 1.0;
wi = 0.0;
for ( m = 1; m < mmax; m += 2 ) {
for ( i = m; i <= n; i += istep )
{
j = i + mmax; tempr = wr * data[ j - 1 ] - wi * data[ j ];
tempi = wr * data[ j ] + wi * data[ j - 1 ];
data[ j - 1 ] = data[ i - 1 ] - tempr;
data[ j ] = data[ i ] - tempi;
data[ i - 1 ] += tempr;
data[ i ] += tempi;
}
wr = ( wtemp = wr ) * wpr - wi * wpi + wr; wi = wi * wpr + wtemp * wpi + wi;
}
mmax = istep;
}
}
The downloaded code was devoid of comments. The author added comments adapted from the ones given in The Reference (p. 513). As indicated in the code, bit-reversal is the process preliminary to the execution of the Danielson-Lanczos algorithm. This process is the main focus of the article.
Instead of running the FFT code on actual data, it is a simple matter to write a tracing program that simulates the bit-reversal process. The following C# console application traces verbatim the logic of the bit-reversal section of the four1
C++ function.
using System;
namespace TraceNRBR
{
class Program
{
static void Main( string[] args )
{
int n, m, j, i;
int nn = 4;
n = nn << 1;
Console.WriteLine( "\nnn = {0}, n == {1}", nn, n );
j = 1;
for ( i = 1; i < n; i += 2 )
{
Console.WriteLine( "\nj == {0}, i == {1}", j, i );
if ( j > i )
{
Console.WriteLine
( "\n\tSWAP( data[ {0} ], data[ {1} ] )", j - 1, i - 1 );
Console.WriteLine( "\tSWAP( data[ {0} ], data[ {1} ] )\n", j, i );
}
m = nn;
while ( m >= 2 && j > m )
{
Console.WriteLine( "m == {0}, j == {1}", m, j );
j -= m;
m >>= 1;
}
j += m;
}
Console.WriteLine();
}
}
}
Upon execution, the program produces the following output:
nn = 4, n == 8
j == 1, i == 1
j == 5, i == 3
SWAP( data[ 4 ], data[ 2 ] )
SWAP( data[ 5 ], data[ 3 ] )
m == 4, j == 5
j == 3, i == 5
j == 7, i == 7
m == 4, j == 7
m == 2, j == 3
Press any key to continue . . .
Considering the four sample complex data points 5 + 4j
, 4 + 3j
, 3 + 2j
, 2 + j
, they would be stored in the data array as follows:
Integer index: 0 1 2 3 4 5 6 7
Binary index: 000 001 010 011 100 101 110 111
Data array: 5 4 4 3 3 2 2 1
R I R I R I R I
From the execution trace, the numbers 3 + 2j and 4 + 3j
would be swapped. However, there are no bit reversals: in the case of integer indices 4
and 2
, 100
is not the reversal of 010
, while in the case of 5
and 3
, 101
is not the reversal of 011
. However, indexing binary pairs of array elements the indices of the four sample complex numbers would be 00
, 01
, 10
, and 11
, respectively. In that case, 01
and 10
indeed are reversals of each other, and the code performs the intended action. Thus the description of bit-reversal at the single array-element indexing level is incorrect. The following figure shows the output from the C# tracing program for 8 hypothetical complex numbers.

Again, the correct indexing of data points is at the complex number level: the complex numbers with indices 001 and 100 (shown in bold on the right-hand side of the figure) would be swapped, as would the numbers with indices 011
and 110
(shown in italics). Clearly the indices in those pairs are bit-reversals of each other. Hence, it would have been better to implement the data array in terms of a struct
or a class
containing the real and imaginary parts of a complex number, and the size of the array would have been the number of complex data points.
C# Implementation of Explicit Bit-Reversal Over Unsigned Integers
Even though their real and imaginary parts can be accessed individually, complex numbers are atomic. As demonstrated by the two executions of the C# tracing program, they are single entities. So, for the purposes of using a C# implementation of explicit bit-reversal, the following minimal class is defined to create and manipulate complex numbers. (The author is aware that C# supports complex numbers, but it does not provide all of the required operations.) The class is implemented as a library, which in turn makes use of some functions in a utility library. The complete code is in the attached ZIP file. All the functionality is properly decribed by comments, and will not be described in detail.
using System;
using System.Text;
using Util_Lib;
namespace ComplexLib
{
public class Complex
{
public static double piDIV2 = Math.PI / (double)2.0;
private static char degChar = '\u00B0';
private double real = 0.0,
imaginary = 0.0;
private double magnitude = 0.0,
phaseDeg = 0.0,
phaseRad = 0.0;
and set its polar form accordingly.
public Complex( double _real, double _imag )
{
real = _real; imaginary = _imag;
magnitude = Math.Sqrt( real * real + imaginary * imaginary );
if ( !UtilFn.NearZero( real ) )
{
phaseRad = Math.Atan( imaginary / real );
}
else
{
phaseRad = imaginary > 0.0 ? piDIV2 : -piDIV2;
}
phaseDeg = UtilFn.RadiansToDegrees( phaseRad );
}
public static Complex PolarToRect( double r, double theta )
{
return new Complex( r * Math.Cos( theta ), r * Math.Sin( theta ) );
}
public Complex Conjugate()
{
return new Complex( real, -imaginary );
}
public override string ToString()
{
return String.Format( "{0,10:000.00000} + ({1,12:0000.00000})i ",
real, imaginary );
}
public string ToLongString( bool newLine = false )
{
StringBuilder sb = new StringBuilder();
sb.Append( ToString() );
if ( newLine )
{
sb.Append( "\n" );
}
sb.Append( String.Format( "== {0,10:000.00000} /_ {1,10:000.00000}{2}\n",
Magnitude, PhaseDegrees, degChar ) );
return sb.ToString();
}
public double Magnitude
{
get { return magnitude; }
}
public double PhaseRadians
{
get { return phaseRad; }
}
public double PhaseDegrees
{
get { return phaseDeg; }
}
public static Complex operator +( Complex a, Complex b )
{
return new Complex( a.real + b.real, a.imaginary + b.imaginary );
}
public static Complex operator -( Complex a, Complex b )
{
return new Complex( a.real - b.real, a.imaginary - b.imaginary );
}
public static Complex operator *( double d, Complex a )
{
return new Complex( d * a.real, d * a.imaginary );
}
public static Complex operator *( Complex a, Complex b )
{
return new Complex( ( a.real * b.real ) - ( a.imaginary * b.imaginary ),
( a.real * b.imaginary ) + ( a.imaginary * b.real ) );
}
public static Complex operator /( Complex a, double d )
{
return new Complex( a.real / d, a.imaginary / d );
}
public static Complex operator /( Complex a, Complex b )
{
Complex bConj = b.Conjugate();
Complex numerator = a * bConj,
denominator = b * bConj;
return numerator / denominator.real;
}
public static Complex[] ToComplexArray( double[] x )
{
int N = x.Length;
Complex[] C = new Complex[ N ];
for ( int i = 0; i < N; ++i )
{
C[ i ] = new Complex( x[ i ], 0.0 );
}
return C;
}
public static Complex Copy( Complex src )
{
return new Complex( src.real, src.imaginary );
}
public static Complex[] Copy( Complex[] src )
{
int n = src.Length;
Complex[] dst = new Complex[ n ];
for ( int i = 0; i < n; ++i )
{
dst[ i ] = Copy( src[ i ] );
}
return dst;
}
public static void Swap( ref Complex x, ref Complex y )
{
Complex t = Copy( x );
x = y;
y = t;
}
}
}
The explicit bit-reversal of an unsigned integer can be implemented by means of shift and bit-manipulation operations, that is, without using division or multiplication operations.
using System;
using System.Text;
using Util_Lib;
namespace ReverseBits
{
static void Main( string[] args )
{
int sigBits;
sigBits = 32;
UInt32 n = 0x700B000A, nRev = RevBits( n );
Display( n, nRev, sigBits );
n = 0xA0B0C0D0; nRev = RevBits( n );
Display( n, nRev, sigBits );
n = 0xB0C0D0E0; nRev = RevBits( n );
Display( n, nRev, sigBits );
n = 0xA1B2C3D4; nRev = RevBits( n );
Display( n, nRev, sigBits );
n = 0xCDEF1235; nRev = RevBits( n );
Display( n, nRev, sigBits );
Process( 3 );
Process( 4 );
Console.WriteLine();
}
private static void Process( int sigBits )
{
List<UInt32> swapped = new List<uint>();
Console.WriteLine( "\n--------------------\nSignificant bits: {0}",
sigBits );
for ( UInt32 n = 0; n < UtilFn.PowerOf2( sigBits ); ++n )
{
UInt32 nRev = RevBits( n );
Display( n, nRev, sigBits, swapped );
}
}
public static UInt32 RevBits( UInt32 n )
{
UInt32 nCopy = n, nRev = 0;
for ( UInt32 i = 0; i < 32; ++i )
{
nRev |= n & mask;
n <<= 1;
nRev >>= 1;
}
nRev <<= 1;
if ( ( nCopy & mask ) >= mask )
{
nRev |= (UInt32)1;
}
return nRev;
}
public static void Display
( UInt32 n, UInt32 nRev, int sigBits, List<UInt32> swapped = null )
{
string nStr = Convert.ToString( n, toBase: 2 ).PadLeft( 32, '0' ),
nRevStr = Convert.ToString( nRev, toBase: 2 ).PadLeft( 32, '0' );
DisplaySplitStr( " n", nStr, n );
DisplaySplitStr( "nRev", nRevStr, nRev );
nRev >>= ( 32 - sigBits );
if ( nRev != n )
{
if ( swapped != null && !swapped.Contains( n ) && !swapped.Contains( nRev ) )
{
Console.WriteLine( "\n\t ==> SWAP( arr[ {0} ], arr[ {1} ]",
n, nRev );
swapped.Add( n );
swapped.Add( nRev );
}
}
Console.WriteLine();
}
public static void DisplaySplitStr( string name, string str, UInt32 number )
{
Console.Write( "\n {0} = 0x{1:X8} = ", name, number );
int j = 0;
for ( int i = 0; i < str.Length; ++i )
{
Console.Write( "{0}", str[ i ] );
++j;
if ( j == 4 )
{
Console.Write( " " ); j = 0;
}
}
}
}
}
Program Output
Upon execution, the following output is produced:
n = 0x700B000A = 0111 0000 0000 1011 0000 0000 0000 1010
n = 0x700B000A = 0111 0000 0000 1011 0000 0000 0000 1010
nRev = 0x5000D00E = 0101 0000 0000 0000 1101 0000 0000 1110
n = 0xA0B0C0D0 = 1010 0000 1011 0000 1100 0000 1101 0000
nRev = 0x0B030D05 = 0000 1011 0000 0011 0000 1101 0000 0101
n = 0xB0C0D0E0 = 1011 0000 1100 0000 1101 0000 1110 0000
nRev = 0x070B030D = 0000 0111 0000 1011 0000 0011 0000 1101
n = 0xA1B2C3D4 = 1010 0001 1011 0010 1100 0011 1101 0100
nRev = 0x2BC34D85 = 0010 1011 1100 0011 0100 1101 1000 0101
n = 0xCDEF1235 = 1100 1101 1110 1111 0001 0010 0011 0101
nRev = 0xAC48F7B3 = 1010 1100 0100 1000 1111 0111 1011 0011
--------------------
Significant bits: 3
n = 0x00000000 = 0000 0000 0000 0000 0000 0000 0000 0000
nRev = 0x00000000 = 0000 0000 0000 0000 0000 0000 0000 0000
n = 0x00000001 = 0000 0000 0000 0000 0000 0000 0000 0001
nRev = 0x80000000 = 1000 0000 0000 0000 0000 0000 0000 0000
==> SWAP( arr[ 1 ], arr[ 4 ]
n = 0x00000002 = 0000 0000 0000 0000 0000 0000 0000 0010
nRev = 0x40000000 = 0100 0000 0000 0000 0000 0000 0000 0000
n = 0x00000003 = 0000 0000 0000 0000 0000 0000 0000 0011
nRev = 0xC0000000 = 1100 0000 0000 0000 0000 0000 0000 0000
==> SWAP( arr[ 3 ], arr[ 6 ]
n = 0x00000004 = 0000 0000 0000 0000 0000 0000 0000 0100
nRev = 0x20000000 = 0010 0000 0000 0000 0000 0000 0000 0000
n = 0x00000005 = 0000 0000 0000 0000 0000 0000 0000 0101
nRev = 0xA0000000 = 1010 0000 0000 0000 0000 0000 0000 0000
n = 0x00000006 = 0000 0000 0000 0000 0000 0000 0000 0110
nRev = 0x60000000 = 0110 0000 0000 0000 0000 0000 0000 0000
n = 0x00000007 = 0000 0000 0000 0000 0000 0000 0000 0111
nRev = 0xE0000000 = 1110 0000 0000 0000 0000 0000 0000 0000
--------------------
Significant bits: 4
n = 0x00000000 = 0000 0000 0000 0000 0000 0000 0000 0000
nRev = 0x00000000 = 0000 0000 0000 0000 0000 0000 0000 0000
n = 0x00000001 = 0000 0000 0000 0000 0000 0000 0000 0001
nRev = 0x80000000 = 1000 0000 0000 0000 0000 0000 0000 0000
==> SWAP( arr[ 1 ], arr[ 8 ]
n = 0x00000002 = 0000 0000 0000 0000 0000 0000 0000 0010
nRev = 0x40000000 = 0100 0000 0000 0000 0000 0000 0000 0000
==> SWAP( arr[ 2 ], arr[ 4 ]
n = 0x00000003 = 0000 0000 0000 0000 0000 0000 0000 0011
nRev = 0xC0000000 = 1100 0000 0000 0000 0000 0000 0000 0000
==> SWAP( arr[ 3 ], arr[ 12 ]
n = 0x00000004 = 0000 0000 0000 0000 0000 0000 0000 0100
nRev = 0x20000000 = 0010 0000 0000 0000 0000 0000 0000 0000
n = 0x00000005 = 0000 0000 0000 0000 0000 0000 0000 0101
nRev = 0xA0000000 = 1010 0000 0000 0000 0000 0000 0000 0000
==> SWAP( arr[ 5 ], arr[ 10 ]
n = 0x00000006 = 0000 0000 0000 0000 0000 0000 0000 0110
nRev = 0x60000000 = 0110 0000 0000 0000 0000 0000 0000 0000
n = 0x00000007 = 0000 0000 0000 0000 0000 0000 0000 0111
nRev = 0xE0000000 = 1110 0000 0000 0000 0000 0000 0000 0000
==> SWAP( arr[ 7 ], arr[ 14 ]
n = 0x00000008 = 0000 0000 0000 0000 0000 0000 0000 1000
nRev = 0x10000000 = 0001 0000 0000 0000 0000 0000 0000 0000
n = 0x00000009 = 0000 0000 0000 0000 0000 0000 0000 1001
nRev = 0x90000000 = 1001 0000 0000 0000 0000 0000 0000 0000
n = 0x0000000A = 0000 0000 0000 0000 0000 0000 0000 1010
nRev = 0x50000000 = 0101 0000 0000 0000 0000 0000 0000 0000
n = 0x0000000B = 0000 0000 0000 0000 0000 0000 0000 1011
nRev = 0xD0000000 = 1101 0000 0000 0000 0000 0000 0000 0000
==> SWAP( arr[ 11 ], arr[ 13 ]
n = 0x0000000C = 0000 0000 0000 0000 0000 0000 0000 1100
nRev = 0x30000000 = 0011 0000 0000 0000 0000 0000 0000 0000
n = 0x0000000D = 0000 0000 0000 0000 0000 0000 0000 1101
nRev = 0xB0000000 = 1011 0000 0000 0000 0000 0000 0000 0000
n = 0x0000000E = 0000 0000 0000 0000 0000 0000 0000 1110
nRev = 0x70000000 = 0111 0000 0000 0000 0000 0000 0000 0000
n = 0x0000000F = 0000 0000 0000 0000 0000 0000 0000 1111
nRev = 0xF0000000 = 1111 0000 0000 0000 0000 0000 0000 0000
Press any key to continue . . .
Conclusion
This article has dealt with the design, implementation and testing of a C# program to explicitly perform the bit-reversal operation required for the computation of the FFT of a set of data. The test runs demonstrated that the program is a correct implementation of bit-reversal over unsigned integer indices.
History
- 28th March, 2023: Initial version