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