This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
FFT differences
- To: gsl-discuss at sourceware dot cygnus dot com
- Subject: FFT differences
- From: Jim Haefner <jhaefner at biology dot usu dot edu>
- Date: Thu, 17 May 2001 13:57:54 -0600
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