Попытка устранить незнакомую ошибку сегментации в C - PullRequest
1 голос
/ 18 февраля 2011

Я пытаюсь узнать, как использовать быстрое преобразование Фурье, и скопировал алгоритм FFT из Числовых реципиентов в C, называемый 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 (законно)

Ответы [ 4 ]

3 голосов
/ 18 февраля 2011

В вашем массиве только 30 значений, но вы говорите four1, что ваш массив имеет длину 64 числа с плавающей запятой.

four1(x-1, n, 1);
1 голос
/ 18 февраля 2011

Вот трассировка стека вашей аварии:

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

Ваш указатель выходного_файла, а также кадр стека были искажены из-за этой строки, которая неверно принимала размер 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/

1 голос
/ 18 февраля 2011

Ваш входной массив с плавающей точкой имеет длину 30, но вы передаете 64 в качестве значения nn.n затем устанавливается в два раза nn, что составляет 128:

n = nn << 1;

(кстати, эта строка является просто запутанным способом записи n = nn * 2, учитывая, что nn это unsigned long).

Затем функция получает доступ к значениям до data[128], что эквивалентно x[127] и далеко за пределами массива.

0 голосов
/ 18 февраля 2011

Будет легче помочь, если вы сообщите нам, какую операционную систему и компилятор вы используете.

При этом ошибка сегментации обычно означает ошибку с указателями.Я предполагаю, что вы скопировали это вручную?Затем найдите крошечную опечатку, что-то вроде *p для **p или что-то в равной степени тривиальное.

Если вы работаете в системе UNIX, задайте для своего corelimit большое значение, запустите его снова.Вы получите coredump, скорее всего.используйте gdb, чтобы выяснить, где происходит настоящая ошибка.

Я бы также проверил строку

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

Обычное соглашение C состоит в том, что индексы находятся в диапазоне от 0 до n-1 , но тест <= предполагает, что я буду в диапазоне до <em>n .

...