Trying to solve an unknwon segmentation fault in C

I'm trying to learn how to use Fast Fourier Transform, and have copied the FFT algorithm from Numerical Recipies in C called four1. I have written a small function test that fourier transforms a simple function. But upon execution the program returns a segmentation fault. I can't find where the fault is can you help me?

#include <stdio.h>
#include <math.h>
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

void four1(float data[], unsigned long nn, int isign)
{
        unsigned long n,mmax,m,j,istep,i;
        double wtemp,wr,wpr,wpi,wi,theta;
        float tempr,tempi;

        n=nn << 1;
        j=1;
        for (i=1;i<n;i+=2) {
                if (j > i) {
                        SWAP(data[j],data[i]);
                        SWAP(data[j+1],data[i+1]);
                }
                m=n >> 1;
                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]-wi*data[j+1];
                                tempi=wr*data[j+1]+wi*data[j];
                                data[j]=data[i]-tempr;
                                data[j+1]=data[i+1]-tempi;
                                data[i] += tempr;
                                data[i+1] += tempi;
                        }
                        wr=(wtemp=wr)*wpr-wi*wpi+wr;
                        wi=wi*wpr+wtemp*wpi+wi;
                }
                mmax=istep;
        }
}

void fourier_transform_test (FILE* output_file)
{
        /* 
        This function serves as a test to see whether my implementation of the
        fft is working or not.
        */
        int n = 30;             // number of samples
        float x[n];             // array that holds all values for x            

        // misc
        int i = 0;

        printf("Running fourier transform t开发者_Go百科ests...\n");
        fprintf(output_file, "# x t\n"); 

        // fill the array x with values to be transformed
        for (i = 0; i <= (n - 1); i++) 
                x[i] = cos((2 * 3.1415 * i) / 10);

        // according to the Numerical Recipies, I have to decrement the pointer to data
        // by one to compensate for the zero-offset
        four1(x-1, 64, 1);

        // loop through the transformed array x and print results to a file
        for (i = 0; i <= (n - 1); i++)
                fprintf(output_file, "%i\t%f\n", i, x[i]);

        fclose(output_file);

}

int main (int argc, char *argv[]) 
{
        // open data_file to write results
        FILE* file; 
        if (argc == 1)
                file = fopen("results.dat", "w");
        else
                file = fopen(argv[1], "w");

        fourier_transform_test(file);

        return 0;
} 

I'm using latest Debian, gcc 4 (pretty sure)

For those of you who want to see the book for yourself: http://www.nrbook.com/a/bookcpdf/c12-2.pdf (its legit)

You only have 30 values in your array but you're telling four1 that your array is 64 floats long.

four1(x-1, n, 1);

Your input array of floats is of length 30, but you pass 64 as the nn value. n is then set to twice nn, which is 128:

n = nn << 1;

(by the way, this line is just an obfuscatory way to write n = nn * 2, given that nn is unsigned long).

The function then accesses values up to data[128], which is equivalent to x[127], and well beyond the bounds of the array.

Here is a stacktrace of your crash:

Breakpoint 1, main (argc=1, argv=0x7fffffffe3c8) at main.c:86
86          if (argc == 1)
(gdb) n
88                  file = fopen("results.dat", "w");
(gdb) n
92          fourier_transform_test(file);
(gdb) print file
$1 = (FILE *) 0x603010
(gdb) n
Running fourier transform tests...

Program received signal SIGSEGV, Segmentation fault.
0x00007ffff7313974 in fclose () from /lib/libc.so.6
(gdb) bt
#0  0x00007ffff7313974 in fclose () from /lib/libc.so.6
#1  0x0000000000400e55 in fourier_transform_test (output_file=0xffffe1e0ffffe1e0) at main.c:78

your output_file pointer and also stack frame was garbaged because of this line, which incorrectly assumed size of x[] to be more than 30 floats:

four1(x-1, 64, 1);

And here is the place which overwrites output_file pointer by reaching over array limits:

Breakpoint 1, main (argc=1, argv=0x7fffffffe3c8) at main.c:90
90          if (argc == 1)
(gdb) n
92                  file = fopen("results.dat", "w");
(gdb) n
96          fourier_transform_test(file);
(gdb) s
fourier_transform_test (output_file=0x603010) at main.c:52
52  {
(gdb) watch output_file 
Hardware watchpoint 2: output_file
(gdb) c
Continuing.
Running fourier transform tests...
Hardware watchpoint 2: output_file

Old value = (FILE *) 0x603010
New value = (FILE *) 0x0
four1 (data=0x7fffffffe1cc, nn=64, isign=1) at main.c:16
16                          SWAP(data[j+1],data[i+1]);
(gdb) bt
#0  four1 (data=0x7fffffffe1cc, nn=64, isign=1) at main.c:16
#1  0x0000000000400dfe in fourier_transform_test (output_file=0x0) at main.c:72
#2  0x0000000000400edb in main (argc=1, argv=0x3f4f07073e9df6ca) at main.c:96
(gdb) l
11          n=nn << 1;
12          j=1;
13          for (i=1;i<n;i+=2) {
14                  if (j > i) {
15                          SWAP(data[j],data[i]);
16                          SWAP(data[j+1],data[i+1]); <---i=39;j=101
17                  }
18                  m=n >> 1;
19                  while (m >= 2 && j > m) {
20                          j -= m;
(gdb) info locals
m = 64
j = 101
istep = 6
wpr = 6.953355807477696e-310
wpi = 0
n = 128
wtemp = 4.9406564584124654e-324
wr = 6.9533490701916841e-310
wi = 6.953355807463467e-310
mmax = 4603839450133099783
i = 39
theta = 6755399441055756
tempr = 0
tempi = -nan(0x7fe260)
(gdb) 

Here is an excellent book about how to implement Fourier Transform - http://www.dspguide.com/

It will be easier to help if you tell us what operating system and compiler you're using.

That said, a segmentation fault usually means an error with pointers. I'm assuming you copied this by hand? Then look for a tiny typo, something like *p for **p or something equally trivial.

If you're on a UNIX system, set your corelimit to a large number, run it again. You'll get a coredump, most likely. use gdb to find out where the actual error occurs.

I'd also check the line

 for (i=m;i<=n;i+=istep) {

Usual C convention is that indices range from 0 to n-1, but a <= test suggests i will range up to n.