尝试解决 C 中未知的分段错误

Question

我正在尝试学习如何使用快速傅立叶变换，并从 C 语言的数值食谱中复制了 FFT 算法，称为 four1。我编写了一个小函数测试，对一个简单函数进行傅里叶变换。但在执行时，程序返回分段错误。我找不到问题出在哪里，你能帮我吗？

#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 tests...\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;
} 

我正在使用最新的 Debian、gcc 4（非常确定）

对于那些想亲自阅读这本书的人： http://www.nrbook.com/a/bookcpdf/c12-2.pdf （其合法性）

Answer 1

您的数组中只有 30 个值，但您告诉 four1 您的数组有 64 个浮点数长。

four1(x-1, n, 1);

Answer 2

您输入的浮点数数组的长度为 30，但您将 64 作为 nn 值传递。然后将 n 设置为两倍 nn，即 128：（

n = nn << 1;

顺便说一下，这一行只是一种混淆的方式来编写 n = nn * 2 code>，假设 nn 是 unsigned long）。

然后，该函数将访问高达 data[128] 的值，这相当于 x[127]，并且远远超出了数组的范围。

Answer 3

这是崩溃的堆栈跟踪：

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

由于这一行，您的output_file指针和堆栈帧被垃圾化，该行错误地假定x[]的大小超过30个浮点：

four1(x-1, 64, 1);

这里是通过到达数组来覆盖output_file指针的地方限制：

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) 

这是一本关于如何实现傅里叶变换的优秀书籍 - http://www.dspguide.com/

Answer 4

如果您告诉我们您正在使用什么操作系统和编译器，我们会更容易为您提供帮助。

也就是说，分段错误通常意味着指针错误。我猜你是手工复制的吧？然后寻找一个微小的拼写错误，例如 **p 的 *p 或同样微不足道的东西。

如果您使用的是 UNIX 系统，请将 corelimit 设置为较大的数字，然后再次运行。您很可能会得到一个核心转储。使用gdb找出实际错误发生的位置。

我还会检查“

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

Usual C convention is thatindexs range from 0 to n-1”这一行，但 <= 测试表明 i 的范围将达到 n。