I implemented a 4-point radix-4 FFT and found that I need to do some manipulation of the output terms to get it to match a dft.

My code is a pretty direct implementation of the matrix-formulation so I'm not clear on what the problem is

//                                |
// radix-4 butterfly matrix form  |  complex multiplication
//                                |
//        +-          -+ +-  -+   |    a+ib
// X[0] = | 1  1  1  1 | |x[0]|   |  * c+id
// X[1] = | 1 -i -1  i | |x[1]|   |    -------
// X[2] = | 1 -1  1 -1 | |x[2]|   |    ac + ibc
// X[3] = | 1  i -1 -i | |x[3]|   |         iad - bd
//        +-          -+ +-  -+   |    ------------------
//                                |

Can anybody spot where I went wrong?

Thanks,

-David

typedef double fp; // base floating-point type

// naiive N-point DFT implementation as reference to check fft implementation against
//
void dft(int inv, struct cfp *x, struct cfp *y, int N) {

long int i, j;
struct cfp w;
fp ang;

for(i=0; i<N; i++) { // do N-point FFT/IFFT
y[i].r = y[i].i = 0;
if (inv) ang =  2*PI*(fp)i/(fp)N;
else     ang = -2*PI*(fp)i/(fp)N;
for (j=0; j<N; j++) {
w.r = cos(j*ang);
w.i = sin(j*ang);
y[i].r += (x[j].r * w.r - x[j].i * w.i);
y[i].i += (x[j].r * w.i + x[j].i * w.r);
}
}

// scale output in the case of an IFFT
if (inv) {
for (i=0; i<N; i++) {
y[i].r = y[i].r/(fp)N;
y[i].i = y[i].i/(fp)N;
}
}

} // dft()

void r4fft4(int inv, int reorder, struct cfp *x, struct cfp *y) {
struct cfp x1[4], w[4];
fp         ang, temp;
int        i;

//                                |
// radix-4 butterfly matrix form  |  complex multiplication
//                                |
//        +-          -+ +-  -+   |    a+ib
// y[0] = | 1  1  1  1 | |x[0]|   |  * c+id
// y[1] = | 1 -i -1  i | |x[1]|   |    -------
// y[2] = | 1 -1  1 -1 | |x[2]|   |    ac + ibc
// y[3] = | 1  i -1 -i | |x[3]|   |         iad - bd
//        +-          -+ +-  -+   |    ------------------
//                                |

if (inv) ang =  2*PI/(fp)4; // invert sign for IFFT
else     ang = -2*PI/(fp)4;
//
w[1].r = cos(ang*1); w[1].i = sin(ang*1); // twiddle1 = exp(-2*pi/4 * 1);
w[2].r = cos(ang*2); w[2].i = sin(ang*2); // twiddle2 = exp(-2*pi/4 * 2);
w[3].r = cos(ang*3); w[3].i = sin(ang*3); // twiddle3 = exp(-2*pi/4 * 3);

//         *1       *1       *1       *1
y[0].r  = x[0].r + x[1].r + x[2].r + x[3].r;
y[0].i  = x[0].i + x[1].i + x[2].i + x[3].i;
//         *1       *-i      *-1      *i
x1[1].r = x[0].r + x[1].i - x[2].r - x[3].i;
x1[1].i = x[0].i - x[1].r - x[2].i + x[3].r;
//         *1       *-1      *1       *-1
x1[2].r = x[0].r - x[1].r + x[2].r - x[3].r;
x1[2].i = x[0].i - x[1].i + x[2].i - x[3].i;
//         *1       *i       *-1      *-i
x1[3].r = x[0].r - x[1].i - x[2].r + x[3].i;
x1[3].i = x[0].i + x[1].r - x[2].i - x[3].r;
//
y[1].r = x1[1].r*w[1].r - x1[1].i*w[1].i; // scale radix-4 output
y[1].i = x1[1].i*w[1].r + x1[1].r*w[1].i;
//
y[2].r = x1[2].r*w[2].r - x1[2].i*w[2].i; // scale radix-4 output
y[2].i = x1[2].i*w[2].r + x1[2].r*w[2].i;
//
y[3].r = x1[3].r*w[3].r - x1[3].i*w[3].i; // scale radix-4 output
y[3].i = x1[3].i*w[3].r + x1[3].r*w[3].i;

// reorder output stage ... mystery as to why I need this
if (reorder) {
temp = y[1].r;
y[1].r = -1*y[1].i;
y[1].i = temp;
//
y[2].r = -1*y[2].r;
//
temp = y[3].r;
y[3].r = y[3].i;
y[3].i = -1*temp;
}

// scale output for inverse FFT
if (inv) {
for (i=0; i<4; i++) { // scale output by 1/N for IFFT
y[i].r = y[i].r/(fp)4;
y[i].i = y[i].i/(fp)4;
}
}

} // r4fft4()
• Can you also show us some sample input and output data for each ? Jul 27 '12 at 15:25
• In addition to the bit-reversal order issue, is there a 2x or 4x difference - some implementations scale the forward fft, some the reverse, and some scale both...
– twalberg
Jul 27 '12 at 15:25
• It's not a re-ordering problem as re-ordering permutes the entries of y as I understand it. I can fix the problem if I change ang = -2*PI; rather than ang = -2*PI/(fp)4; I don't need to re-order the terms and my consistency test versus the dft passes with 0 errors. I think this is equivalent to a 90deg phase shift for the twiddle factors. However, this doesn't seem consistent with the mathematics ... what am i missing?
– user1211582
Jul 27 '12 at 15:41

I've just ported a radix-4 DIF fft from S. Burrus Fortran code to Java. Actually it lacks several optimization, first of all the table driven twiddle factor (sin and cos factors should be pre-calculated). This should speed the fft a bit more (maybe 50%). I've to hack a bit for that but if someone has the correct answer I'll be very happy and grateful. I'll post the optimized code asap I hope maybe with some speed tests vs radix-2 algorithm.

More, the multiplications by 1 and sqrt(-1) aren't removed. Removing them will speed a little more. But overall IMHO radix-4 seems be not more than 25% faster than a radix-2, so I don't know if the speed/complexity ratio is really worth. Keep in mind that very optimized libraries like FFTW are largely available and used, so this effort could be just a personal 'divertissment'!

Here is the java code. Porting it to C, C++ or C# should be very easy.

public static void FFTR4(double[] X, double[] Y, int N, int M) {
// N = 4 ^ M
int N1,N2;
int I1, I2, I3;
double CO1,CO2,CO3,SI1,SI2,SI3;
double A,B,C,E;
double R1,R2,R3,R4;
double S1,S2,S3,S4;
// N = 1 << (M+M);
N2 = N;
I2 = 0; I3 = 0;
for (int K=0; K<M; ++K) {
N1 = N2;
N2 = N2 / 4;
E = PI2 / (double)N1;
A = 0.0;
for (int J=0; J < N2; ++J) {
A = J*E;
B = A + A;
C = A + B;
//Should be pre-calculated for optimization
CO1 = Math.cos(A);
CO2 = Math.cos(B);
CO3 = Math.cos(C);
SI1 = Math.sin(A);
SI2 = Math.sin(B);
SI3 = Math.sin(C);
for (int I = J; I<N; I+=N1) {
I1 = I + N2;
I2 = I1 + N2;
I3 = I2 + N2;
R1 = X[I] + X[I2];
R3 = X[I] - X[I2];
S1 = Y[I] + Y[I2];
S3 = Y[I] - Y[I2];
R2 = X[I1] + X[I3];
R4 = X[I1] - X[I3];
S2 = Y[I1] + Y[I3];
S4 = Y[I1] - Y[I3];
X[I] = R1 + R2;
R2 = R1 - R2;
R1 = R3 - S4;
R3 = R3 + S4;
Y[I] = S1 + S2;
S2 = S1 - S2;
S1 = S3 + R4;
S3 = S3 - R4;
X[I1] = CO1*R3 + SI1*S3;
Y[I1] = CO1*S3 - SI1*R3;
X[I2] = CO2*R2 + SI2*S2;
Y[I2] = CO2*S2 - SI2*R2;
X[I3] = CO3*R1 + SI3*S1;
Y[I3] = CO3*S1 - SI3*R1;
}
}
}

double T;
int J = 0;
N2 = N>>2;
for (int I=0; I < N-1; I++) {
if (I < J) {
T = X[I];
X[I] = X[J];
X[J] = T;
T = Y[I];
Y[I] = Y[J];
Y[J] = T;
}
N1 = N2;
while ( J >= 3*N1 ) {
J -= 3*N1;
N1 >>= 2;
}
J += N1;
}
}

Here is the original Radix-4 DIF FORTRAN code by Sidney Burrus:

First, your supposed 'radix-4 butterfly' is a 4 point DFT, not an FFT. It has 16 complex (ie: N squared) operations. A typical 4 point FFT would have only Nlog(base 2)N (= 8 for N = 4). Second, you have some supposed w[ ].r and w[ ].i 'scale' factors that don't belong. Perhaps you obtained them from a radix-4 butterfly shown in a larger graph. Such a butterfly would have some interstage twiddles appended to it, but they are not actually part of the butterfly. A 4 point FFT only has an internal butterfly of -j when designed for a negative exponent FFT.

Rather than try to fix your code, it's just as easy to write my own, as shown below (DevC++ compiler; outputs appended at end of code):

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <cmath>
using namespace std;
void fft4(double* r, double* i);    // prototype declaration
int main (int nNumberofArgs, char* pszArgs[ ] ) { // arguments needed for Dev C++ I/O

double r[4] = {1.5, -2.3, 4.65, -3.51}, i[4] = {-1.0, 2.6, 3.75, -2.32} ;
long n, k, j;      double  yr[4] = {0.}, yi[4] = {0.};
double ang, C, S, twopi = 6.2831853071795865;

cout<<"\n original real/imag data";
cout<<"\n n         r[n]            i[n]\n";
for (n = 0; n < 4; n++)  {
printf("%2d\t%9.4f\t%9.4f\n",n,r[n],i[n]);
} //end for loop over n

// 4 point DFT
for (k = 0; k < 4; k++) {
ang = twopi*k/4;
for (j = 0; j < 4; j++) {
C = cos(j*ang);       S = sin(j*ang);
yr[k] = yr[k] + r[j]*C + i[j]*S;   // ( C - jS )*( r + ji )
yi[k] = yi[k] + i[j]*C - r[j]*S;   // = ( rC + iS ) + j( iC - rS )
}
}

cout<<"\n 4 point DFT results";
cout<<"\n n         yr[n]           yi[n]           amplitude       phase(radians)\n";
double amp, phase;
for (n = 0; n < 4; n++)  {
yr[n] = yr[n]/4 ;      yi[n] = yi[n]/4 ;  // scale outputs
amp = sqrt( yr[n]*yr[n] + yi[n]*yi[n] ) ;
phase = atan2( yi[n], yr[n] ) ;
printf("%2d\t%9.4f\t%9.4f\t%9.4f\t%9.4f\n",n,yr[n],yi[n],amp,phase);
} //end for loop over n

fft4(r, i) ;

cout<<"\n 4 point FFT results";
cout<<"\n n         r[n]            i[n]            amplitude       phase(radians)\n";

for (n = 0; n < 4; n++)  {
r[n] = r[n]/4 ;      i[n] = i[n]/4 ;  // scale outputs
amp = sqrt( r[n]*r[n] + i[n]*i[n] ) ;
phase = atan2( i[n], r[n] ) ;
printf("%2d\t%9.4f\t%9.4f\t%9.4f\t%9.4f\n",n,r[n],i[n],amp,phase);
} //end for loop over n

fft4(i, r); // this is an inverse FFT (complex in/out routine)

cout<<"\n 4 point inverse FFT results";
cout<<"\n n         r[n]            i[n]\n";
for (n = 0; n < 4; n++)  {
printf("%2d\t%9.4f\t%9.4f\n",n,r[n],i[n]);
} //end for loop over n

system ("PAUSE");
return 0;
} // end main
//************************ fft4 **********
void fft4(double* r, double* i) {
double t;

t = r[0]; r[0] = t + r[2]; r[2] = t - r[2];
t = i[0]; i[0] = t + i[2]; i[2] = t - i[2];
t = r[1]; r[1] = t + r[3]; r[3] = t - r[3];
t = i[1]; i[1] = t + i[3]; i[3] = t - i[3];

t = r[3]; r[3] = i[3]; i[3] = -t; // (r + ji)*(-j)

t = r[0]; r[0] = t + r[1]; r[1] = t - r[1];
t = i[0]; i[0] = t + i[1]; i[1] = t - i[1];
t = r[2]; r[2] = t + r[3]; r[3] = t - r[3];
t = i[2]; i[2] = t + i[3]; i[3] = t - i[3];

t = r[1]; r[1] = r[2]; r[2] = t;  // swap 1
t = i[1]; i[1] = i[2]; i[2] = t;  //  and 2
} // end fft4

original real/imag data
n         r[n]            i[n]
0         1.5000         -1.0000
1        -2.3000          2.6000
2         4.6500          3.7500
3        -3.5100         -2.3200

4 point DFT results
0         0.0850          0.7575          0.7623          1.4591
1         0.4425         -1.4900          1.5543         -1.2821
2         2.9900          0.6175          3.0531          0.2037
3        -2.0175         -0.8850          2.2031         -2.7282

4 point FFT results
0         0.0850          0.7575          0.7623          1.4591
1         0.4425         -1.4900          1.5543         -1.2821
2         2.9900          0.6175          3.0531          0.2037
3        -2.0175         -0.8850          2.2031         -2.7282

4 point inverse FFT results
n         r[n]            i[n]
0         1.5000         -1.0000
1        -2.3000          2.6000
2         4.6500          3.7500
3        -3.5100         -2.3200

First, the input data (4 real, 4 imaginary) is printed out. Then a 4 point DFT is taken. The results (yr[ ] and yi[ ] plus amp/phase) are printed out. Since the r[ ] and i[ ] original data were not overwritten when doing the DFT, those inputs are reused as inputs to the 4 point FFT. Note that the latter has fewer +/- operations than the DFT.

The code for the FFT is not particularly elegant nor efficient – there are many ways of doing butterflies. The code above corresponds to the four radix-2 butterflies shown in Rabiner and Gold's book “Theory and Application of Digital Signal Processing” (p. 580, Fig. 10.9), with twiddles modified to reflect a negative exponent (the ones used for the figure in the book were positive). Note that there is only one twiddle of -j in the code, and this does not require a multiply (it's a swap/sign change).

After the FFT, the results are printed. They're the same as the DFT

And finally, the scaled results from the FFT are used as inputs to an inverse FFT. This is accomplished via the 'exchange' or 'reverse the list' method (ie: if FFT(r,i) is a forward FFT, then FFT(i,r) is an inverse – provided, of course, that the FFT is capable of handling complex inputs/outputs – in other words – no 'real only' routines, which usually presume that imaginary inputs are zero). This method was described nearly 25 years ago in:

P. Duhamel, B. Piron, J. M. Etcheto, “On Computing the Inverse DFT,” IEEE Transactions on Acoustics, Speech and Signal Processing, vol. 36, Feb. 1988, pp. 285-286.

The result of the inverse is then printed out. It is the same as the original input data.