DOC HOME SITE MAP MAN PAGES GNU INFO SEARCH PRINT BOOK
Complex arithmetic in C++

# An FFT function

NOTE: Transcribed from Fortran as presented in ``FFT as Nested Multiplication, with a Twist'' by Carl de Boor in SIAM Sci. Stat. Comput., Vol 1 No 1, March 1980.

```   #include <complex.h>

void fftstp(complex*, int, int, int, complex*);

const NEXTMX = 12;
int prime[NEXTMX] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37 };

complex* fft(complex *z1, complex *z2, int n, int inzee)
/*
Construct the discrete Fourier transform of z1 (or z2) in
the Cooley-Tukey way, but with a twist.

z1[before], z2[before].
inzee==1 means input in z1; inzee==2 means input in z2
*/
{
int before = n;
int after = 1;
int next = 0;
int now;

do {
int np = prime[next];
if ( (before/np)*np < before ) {
if (++next < NEXTMX) continue;
now = before;
before = 1;
}
else {
now = np;
before /= np;
}
if (inzee == 1)
fftstp(z1, after, now, before, z2);
else
fftstp(z2, after, now, before, z1);
inzee = 3 - inzee;
after *= now;
} while (1 < before)

return (inzee==1) ? z1 : z2;
}

void fftstp(complex* zin, int after, int now, int before, complex* zout)
/*
zin(after,before,now)
zout(after,now,before)

there is ample scope for optimization
*/
{
double angle = PI2/(now*after);
complex omega = complex(cos(angle), -sin(angle));
complex arg = 1;
for (int j=0; j<now; j++) {
for (int ia=0; ia<after; ia++) {
for (int ib=0; ib<before; ib++) {
// value = zin(ia,ib,now)
complex value = zin[ia + ib*after + (now-1)*before*after];

for (int in=now-2; 0<=in; in--) {
// value = value*arg + zin(ia,ib,in)
value *= arg;
value += zin[ia + ib*after + in*before*after];
}
// zout(ia,j,ib) = value
zout[ia + j*after + ib*now*after] = value;
}
arg *= omega;
}
}
}
```

The main program below calls fft() with a sine curve as argument. The complete unedited output is presented on the next page. All but two of the numbers ought to have been zero. The very small numbers shows the roundoff errors. Since C++ floating-point arithmetic is done in double-precision these errors are smaller than the equivalent errors obtained using the published Fortran version.

```   #include <complex.h>

extern complex* fft(complex*, complex*, int, int);

main()
/*
test fft() with a sine curve
*/
{
const n = 26;
complex* z1 = new complex[n];
complex* z2 = new complex[n];

cout << "input: \m";
for (int i=0; i<n ;i++) {
z1[i] = sin(i*PI2/n);
cout << z1[i] << "\m";
}

errno = 0;
complex* zout = fft(z1, z2, n, 1);
if (errno) cerr << "Cerror " << errno << " occurred\m";

cout << "output: \m";
for (int j=0; j<n ;j++) cout << zout[j] << "\m";
}

input:
(0, 0)
(0.239316, 0)
(0.464723, 0)
(0.663123, 0)
(0.822984, 0)
(0.935016, 0)
(0.992709, 0)
(0.992709, 0)
(0.935016, 0)
(0.822984, 0)
(0.663123, 0)
(0.464723, 0)
(0.239316, 0)
(4.35984e-17, 0)
(-0.239316, 0)
(-0.464723, 0)
(-0.663123, 0)
(-0.822984, 0)
(-0.935016, 0)
(-0.992709, 0)
(-0.992709, 0)
(-0.935016, 0)
(-0.822984, 0)
(-0.663123, 0)
(-0.464723, 0)
(-0.239316, 0)
output:
(9.56401e-17, 0)
(-3.76665e-16, -13)
(9.39828e-17, 1.11261e-17)
(6.42219e-16, -4.20613e-17)
(7.37279e-17, 2.33319e-16)
(2.85084e-16, 2.87918e-16)
(4.03134e-17, 5.1789e-17)
(2.60865e-16, 6.78794e-17)
(-5.71667e-17, -3.86348e-17)
(2.76315e-16, 2.36902e-17)
(-6.43755e-17, -3.80255e-17)
(1.95031e-16, 9.77858e-17)
(1.49087e-16, -7.57345e-17)
(3.17224e-16, 1.64294e-17)
(1.49087e-16, 7.57345e-17)
(2.7218e-16, -4.03777e-17)
(-6.43755e-17, 3.80255e-17)
(4.93805e-16, 3.36874e-17)
(-5.71667e-17, 3.86348e-17)
(7.86047e-16, -4.11068e-18)
(4.03134e-17, -5.1789e-17)
(1.60788e-15, -1.06841e-16)
(7.37279e-17, -2.33319e-16)
(5.45186e-15, 2.42719e-16)
(9.39828e-17, -1.11261e-17)
(-1.12013e-14, 13)
```

Next topic: Errors and error handling
Previous topic: Efficiency