FFT differences

[Jim Haefner]

I am just starting to use gsl and see discrepancies between the results of
gsl, numerical recipes,and matlab. 
I realize there will be differences due to algorithms, but it seems to me
that (a) the differences are large,
and (b) the fact that numerical recipes and matlab agree casts aspersions
on gsl.  I cut and pasted the
example code from gsl-ref.info for FFT complex example program.  This code
does a forward FFT on a 128 vector
with real values of 1.0 in the first 10 positions and in the last 10
positions.  Here are the first 10 sets of
coefficients from gsl, num. recipes, and matlab;   I used gsl version 0.7
from the redhat rpm distribution,
matlab 5.3.0.620a (R11), and an unknown version of num. recipes.  Below is
the code used for the gsl results;
below that is the make file.  Numerical recipes and gsl versions were
compiled with gcc version 2.96 in Linux
2.2x from Redhat7.0. 

Any and all corrections to my use of gsl would be greatly appreciated.  At
the moment, I am inclined to
believe Numerical Recipes over gsl, although I would prefer to use gsl.
Does any one know something about
Num.Rec. and Matlab that suggests their code is wrong?

GSL:  Unscaled: (don't divide by sqrt(128))
2.100000e+01 0.000000e+00
2.008450e+01 1.071924e-15
1.748052e+01 1.343489e-15
1.358941e+01 1.525445e-15
8.997623e+00 1.591609e-16
4.370514e+00 -5.305814e-17
3.344068e-01 2.411266e-16
-2.629892e+00 1.644125e-16
-4.261973e+00 -5.905107e-16
-4.551734e+00 -8.908686e-16

NUMERICAL RECIPIES:  Unscaled: 
 2.000000e+01 4.713967e-01
 1.920257e+01 8.314696e-01
 1.692495e+01 9.951847e-01
 1.349139e+01 9.238795e-01
 9.380306e+00 6.343933e-01
 5.143525e+00 1.950903e-01
 1.315192e+00 -2.902847e-01
 -1.672952e+00 -7.071068e-01
 -3.554866e+00 -9.569403e-01
 -4.261449e+00 -9.807853e-01

MATLAB:
20.00000000000000                     
19.20257417032812 + 0.47139673682600i
16.92495277438511 + 0.83146961230254i 
13.49139017453587 + 0.99518472667220i
 9.38030631121151 + 0.92387953251129i  
 5.14352486328997 + 0.63439328416364i
 1.31519212062881 + 0.19509032201613i 
-1.67295181692755 - 0.29028467725446i
-3.55486584620912 - 0.70710678118655i 
-4.26144883171454 - 0.95694033573221i
 



==================
GSL code
==================
/* double version */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_complex.h>
   
#define REAL(z,i) ((z)[2*(i)])
#define IMAG(z,i) ((z)[2*(i)+1])

#define NUMPT (unsigned long)128

int main()
{
    int i;
    /* example from gsl info */
   double data[2*128];
   for (i = 0; i < 128; i++)
     {
        REAL(data,i) = 0.0;
        IMAG(data,i) = 0.0;
     }
   REAL(data,0) = 1.0;
   for (i = 1; i <= 10; i++)
     {
        REAL(data,i) = REAL(data,128-i) = 1.0;
     }
   for (i = 0; i < 128; i++)
     {
       printf ("%d %e %e\n", i, REAL(data,i), IMAG(data,i));
     }
   printf ("\n");
   gsl_fft_complex_radix2_forward (data, 1, 128);
   printf ("Unscaled: \n");
   for (i = 0; i < 128; i++)
     {
       printf ("%d %e %e\n", i, REAL(data,i),
                                IMAG(data,i));
     }
   printf ("\nScaled\n");
   for (i = 0; i < 128; i++)
     {
       printf ("%d %e %e\n", i, REAL(data,i)/sqrt(128),
                                IMAG(data,i)/sqrt(128));
     }
    return 0;
}

=======================
GSL make
=======================

# Gnu C make fft

CFLAGS=-g -Wall
###CFLAGS= -Wall -O2

COMPILER=gcc

##.c.o:
##    $(COMPILER) $(CFLAGS) -c $<

gsl_tfftd:     gsl_tfftd.o   
	$(COMPILER) $(CFLAGS) -o gsl_tfftd  gsl_tfftd.o  -lm -lgsl
-lgslblas

gsl_tfftd.o:   gsl_tfftd.c 
	$(COMPILER) $(CFLAGS) -c gsl_tfftd.c







-- 
James W. Haefner            
Department of Biology   Email: jhaefner@biology.usu.edu
Utah State University   Voice: 435-797-3553
Logan, UT 84322-5305      Fax: 435-797-1575