Working up to an FFT

This is a discussion on Working up to an FFT within the C Programming forums, part of the General Programming Boards category; Hi everyone,
The community here has been very helpful to me with my previous problems, so I hope that you ...

Working up to an FFT

Hi everyone,

The community here has been very helpful to me with my previous problems, so I hope that you can offer some advice with this.

I've been set a problem involving a Fast Fourier Transform (FFT). At the moment I'm simply trying to get the data read in from the data file and output to another file, but I'm getting a segmentation fault when I try to read the data. Here's the code:

The red portion of code may seem odd, but it's a precursor to the FFT. I need an array twice the size of the number of real points. However I've been told to carry out all computations in double precision, but the data in the file is in float. I'm not really sure how I to cast in this case using pointers. (Blue section)

I'm printing out all the data that is read in so that I can compare it to the output file at the end. But I'm getting a segmentation fault partway through. Any pointers (no pun intended) about what I'm doing wrong?

well a couple of weeks have passed, and due to a group project that I'm working on for my course the programming has taken a bit of a back seat. I've not really progressed, and I'm still having problems.

The first problem involves the array sizes. I'm reading the data into an array, and then I have to increase the size of the array to an integer poer of 2 (2^x) so that I can fourier transform it. But I'm having difficulties with the passing of the array size by reference (highlighted in red). Thedata is read into an arrayof size 3792, but the array size taken by padme is 137366! I know it's a problem with my use of pointers, but I can't see what it is and 've been trying for several days now. I'm also getting a gcc compiler error: pointer targets in passing arg 2 of 'padme' differ in signedness.

Code:

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include <sys/types.h>
#include <sys/stat.h>
/*#include "fft.c"*/
double* padme(double *ptr, unsigned *long two_numdata)
{
unsigned long numpads=0, padsize, arraysize;
unsigned int i;
double *fakeptr;
*two_numdata=arraysize;
/*check to see whether the number of data is an integer power of two*/
if(((arraysize-1) & arraysize) == 0)
{
printf("The array size is %ld\n", arraysize);
printf("This is an integer power of 2.\n");
printf("No padding was carried out.\n");
}
else
{
printf("The array has size %ld\n", arraysize);
printf("This is not an integer power of 2, so padding will be carried out\n\n");
padsize=arraysize+numpads; /*increase the size of the array until it
is an integer poer of two by incrementing
the number of padding points that are
required, and testing the total number of points.*/
while(((padsize-1) & padsize) != 0)
{
numpads++;
padsize=arraysize+numpads;
}
printf("final array size = %ld\n", padsize);
/*Reallocate memory to the new array based on the new size*/
ptr=(double*)realloc(ptr, padsize*sizeof(double));
fakeptr=ptr+arraysize;
/*Pad data with 0s to the power of 2*/
for(i=arraysize+1; i<=padsize; i++)
{
*fakeptr=0.0;
fakeptr++;
}
}
return ptr;
}
/*Code for main function adapted from: px390_prj_wrapper.c, M.G.Dowsett 2007, with permission*/
int main()
{
off_t filesize, numdata, two_numdata; /* off_t type definition is
~ unsigned long int and is in sys/types.h*/
FILE *fp, *sp; /* pointer to a file */
float temp;
double *spectrumptr; /* define a pointer to a double type to store the data in */
double *outputptr; /*define a pointer for the output array */
double *dummyptr; /* a spare pointer for use in loops */
int i, isign;
char confirm;
char datafile[]="prj4data.dat";
char *df=&datafile[0]; /* df points to the file name*/
char outfile[]="prj4fft.txt";
/* define a structure to read the file attributes - type definition stat is in sys/stat.h */
/* see http://www.gnu.org/software/libc/manual/ */
struct stat st;
struct stat *buf;
buf=&st;
/* attempt to read the attributes and allocate space for the data */
if (stat(df,buf)==0)
{
filesize=st.st_size; /*size in bytes in bytes*/
numdata=filesize/4; /*number of real floats to read*/
printf("Number of real data is %ld\n", numdata);
two_numdata=2*numdata;
spectrumptr=(double*) malloc(two_numdata*sizeof(double)); /*request a pointer to enough
bytes for the complex array */
if(spectrumptr==NULL)
{
printf("Error allocating memory. The program will now terminate");
exit(EXIT_FAILURE);
}
dummyptr=spectrumptr;
if ((fp=fopen(datafile,"rb"))!=NULL) /*allocate a pointer to the file
and if all is OK read it*/
{
for (i=1; i<=two_numdata; i++)
{
if((i & 1) ==0) /*Check to see if current array element is even or odd.
if even, element is imaginary and is set to 0*/
{
temp=0.0;
*dummyptr=temp;
}
else
{
fread(&temp,sizeof(float),1, fp);
*dummyptr=temp;
}
dummyptr++;
}
printf("Data read into complex array of size %ld.\n", two_numdata);
printf("Imaginary elements have been set to 0.0\n\n");
/*I have added an extra blank line to make the output easier to read and follow*/
fclose(fp);
}
else
{
printf("Cannot open data file. \n");
exit(0);
}
}
else
{
printf("Cannot read file attributes. \n");
exit(0);
}
/*call padme*/
printf("The size of the input array will now be altered if necessary\n");
spectrumptr=padme(spectrumptr, &two_numdata);
/*Call padme again for the output array*/
printf("The size of the output array will now be latered if necessary\n");
outputptr=padme(outputptr, &two_numdata);
/*Nested loops to determine whether a forward or inverse transform is to be carried out.
The loops use Y/N choices to allow the user to decide what they wish to do*/
printf("Do you wish to carry out a Fourier transform, y/n?\n");
confirm=toupper(getchar()); /*Converts the input to upper case to catch, for example, both y and Y*/
while(isign!=1 && isign!=0 && isign!=-1)
{
if(confirm=='Y')
isign=1;
else if(confirm=='N')
{
/*If a forward transform is not required, check to see whether an inverse transform is required*/
printf("Do you wish to carry out an inverse fourier transform, y/n?\n");
confirm=toupper(getchar());
while(isign!=0 && isign!=-1)
{
if(confirm=='Y')
isign=-1;
else if(confirm =='N')
isign=0;
else
printf("Do you wish to carry out an inverse Fourier transform, y/n?");
}
}
else
printf("Do you wish to carry out a Fourier transform, y/n?");
}
printf("\n"); /*Puts an empty line between previous output and following output.
Makes final output easier to read and follow*/
if(isign==1)
{
printf("A Fourier transform of the data will now be carried out\n");
/*mycpxfft(spectrumptr, outputptr, numtotaldata, isign);*/
}
else if(isign==-1)
{
printf("An inverse Fourier transform will now be carried out\n");
/*mycpxfft(spectrumptr, outputptr, numtotaldata, isign);*/
}
else if(isign==0) /*If neither transform is required, the program will simply output the padded input array*/
printf("No transforming of the data will be carried out\n");
dummyptr=spectrumptr;
if ((sp=fopen(outfile,"w")) !=NULL)
{
for(i=1; i<=two_numdata; i++)
{
fprintf(sp,"%f \n",*dummyptr); /* NOTE: \n only puts in a line feed (chr = 10) not the
carriage-return linefeed pair (13 10) required by windows.
Excel/origin will still read*/
dummyptr++;
}
fclose(sp);
}
else
{
printf("Cannot write the text output file. \n");
exit(EXIT_FAILURE);
}
printf("Task done - all OK \n");
free(spectrumptr);
free(outputptr);
return(0);
}

The second problem deals with the portion of code in blue. No matter what I enter for the check (y, n, p etc) it just seems to skip the loops and set isign to 0. the compiler error I get is "implicit declaration of toupper", but I've used this before and it worked fine.

Ok, I've had a breakthrough - strange how it happens isn't it. I've now fixed the problem with the array sizes. I realised that
*two_numdata=arraysize;
was the wrong way round for what I wanted. I also added an extra line at the end of padme
*two_numdata=padsize;
so that when I printout the array at the end of main it prints the right number of data.

The second parameter to padme does have to be a pointer, as that's what my specification says and I'm not allowed to alter it. Also, if I try your suggestions I no longer obtain the correct answers for the padding.

I hadn't noticed the non-initialisation of isign, and didn't realise that toupper was in <ctype.h> - the book I have told me it was in stdlib, but oh well. I'm no longer getting error messages for that. As far as I know (I checked the function listing on this website), getchar works with chars too, so I've left it as is since I find it easier to follow.

I've now added in the fft function, but I'm getting more problems.
1: I'm still suffering from the warning about the signedness of the padme argument.
2: When the program runs the section in red, if I enter 'N' for the fourier transform question it prints the next statement twice. The same thing happens if I input something that's not 'N' or 'Y'. Any ideas why? I've tried chaging the loops several times, and I can't see why it would do so.
3: As soon as the program prints the statement saying what form of fourier transform will be carried out, I get a segmentation fault. No idea why.

Code:

/*
This program reads in the real data from a datafile, and allocates it to every other element in an array of size
2*number of data. The program then checks the size of the array and pads it with additional 0.0s until it is an
integer power of 2 using the function padme. The program then checks to see whether a forward or reverse fourier
transform is required, and carries out the appropriate transformation using the function mycpfft. Finally, the data in
the output array is printed into a text file for processing.
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <ctype.h>
double* padme(double *ptr, unsigned long *two_numdata)
{
unsigned long numpads=0, padsize, arraysize;
unsigned int i;
double *fakeptr;
arraysize=*two_numdata;
/*check to see whether the number of data is an integer power of two*/
if(((arraysize-1) & arraysize) == 0)
{
printf("The array size is %ld\n", arraysize);
printf("This is an integer power of 2.\n");
printf("No padding was carried out.\n\n");
}
else
{
printf("The array size is %ld\n", arraysize);
printf("This is not an integer power of 2, so padding will be carried out\n");
padsize=arraysize+numpads; /*increase the size of the array until it
is an integer poer of two by incrementing
the number of padding points that are
required, and testing the total number of points.*/
while(((padsize-1) & padsize) != 0)
{
numpads++;
padsize=arraysize+numpads;
}
printf("Final array size = %ld\n\n", padsize);
/*Reallocate memory to the new array based on the new size*/
ptr=(double*)realloc(ptr, padsize*sizeof(double));
fakeptr=ptr+arraysize;
/*Pad data with 0s to the power of 2*/
for(i=arraysize+1; i<=padsize; i++)
{
*fakeptr=0.0;
fakeptr++;
}
}
*two_numdata=padsize;
return ptr;
}
/*Function for computing the fast Fourier transform is adapted from: W.H.Press, S.A.Teukolsky, W.T.Vetterling and B.P.Flannery,
Numerical Recipes in C: 2nd edition, Cambridge University Press(1992)*/
void mycpxfft(double data[ ], double fftdata[ ], unsigned long numcomplex, int isign)
{
unsigned long mmax, m, j, istep, i, N;
double wtemp, wr, wpr, wpi, wi, theta;
/*Double precision for the trigonometric recurrences*/
double tempreal, tempim;
N=numcomplex<<1;
j=1;
for(i=1;i<N;i+=2)
{ /*This is the bit-reversal section of the routine*/
if (j>i)
{
/*Exchange the two complex numbers.*/
fftdata[j]=data[i];
fftdata[i]=data[j];
fftdata[j+1]=data[i+1];
fftdata[i+1]=data[j+1];
}
m=numcomplex;
while(m>=2 && j>m)
{
j-= m;
m >>= 1;
}
j+= m;
}
/*Here begins the Danielson-Lanczos section of the routine.*/
mmax=2;
while(N > mmax)
{
/*Outer loop executed log2*nuncomplex times.*/
istep=mmax<<1;
theta=isign*(6.28318530717959/mmax); /*Initialize the trigonometric recurrence.*/
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; /*This is the Danielson-Lanczos formula*/
tempreal=wr*fftdata[j]-wi*fftdata[j+1];
tempim=wr*fftdata[j+1]+wi*fftdata[j];
fftdata[j]=fftdata[i]-tempreal;
fftdata[j+1]=fftdata[i+1]-tempim;
fftdata[i] += tempreal;
fftdata[i+1] += tempim;
}
wr=(wtemp=wr)*wpr-wi*wpi+wr; /*Trigonometric recurrence.*/
wi=wi*wpr+wtemp*wpi+wi;
}
mmax=istep;
}
}
/*Code for main function adapted from: px390_prj_wrapper.c, M.G.Dowsett 2007, with permission*/
int main()
{
off_t filesize, numdata, two_numdata, numcomplex; /* off_t type definition is ~ unsigned long int and is
in sys/types.h*/
FILE *fp, *mp, *sp; /* pointer to a file */
float temp;
double *spectrumptr; /* define a pointer to a double type to store the data in */
double *outputptr; /*define a pointer for the output array */
double *dummyptr; /* a spare pointer for use in loops */
int i, isign=2;
char confirm;
char datafile[]="prj4data.dat";
char *df=&datafile[0]; /* df points to the file name*/
char midfile[]="prj4input.txt"; /*file into which the data can be read after padding,
in order to facilitate checking*/
char outfile[]="prj4fft.txt";
/* define a structure to read the file attributes - type definition stat is in sys/stat.h */
/* see http://www.gnu.org/software/libc/manual/ */
struct stat st;
struct stat *buf;
buf=&st;
/* attempt to read the attributes and allocate space for the data */
if (stat(df,buf)==0)
{
filesize=st.st_size; /*size in bytes in bytes*/
numdata=filesize/4; /*number of real floats to read*/
printf("Number of real data is %ld\n", numdata);
two_numdata=2*numdata;
spectrumptr=(double*) malloc(two_numdata*sizeof(double)); /*request a pointer to enough
bytes for the complex array */
if(spectrumptr==NULL)
{
printf("Error allocating memory. The program will now terminate");
exit(EXIT_FAILURE);
}
dummyptr=spectrumptr;
if ((fp=fopen(datafile,"rb"))!=NULL) /*allocate a pointer to the file
and if all is OK read it*/
{
for (i=1; i<=two_numdata; i++)
{
if((i & 1) ==0) /*Check to see if current array element is even or odd.
if even, element is imaginary and is set to 0*/
{
temp=0.0;
*dummyptr=temp;
}
else
{
fread(&temp,sizeof(float),1, fp);
*dummyptr=temp;
}
dummyptr++;
}
printf("Data read into complex array of size %ld.\n", two_numdata);
printf("Imaginary elements have been set to 0.0\n\n");
/*I have added an extra blank line to make the output easier to read and follow*/
fclose(fp);
}
else
{
printf("Cannot open data file. \n");
exit(0);
}
}
else
{
printf("Cannot read file attributes. \n");
exit(0);
}
/*call padme*/
printf("The size of the input array will now be altered if necessary\n");
spectrumptr=padme(spectrumptr, &two_numdata);
dummyptr=spectrumptr;
printf("The padded input data will now be transferred to a file to preserve it, and allow for error checking.\n");
/*Read the padded data into a text file so that the transformed data can be checked against it*/
if ((mp=fopen(midfile,"w")) !=NULL)
{
for(i=1; i<=two_numdata; i++)
{
fprintf(mp,"%f \n",*dummyptr); /* NOTE: \n only puts in a line feed (chr = 10) not the
carriage-return linefeed pair (13 10) required by windows.
Excel/origin will still read*/
dummyptr++;
}
fclose(mp);
}
else
{
printf("Cannot write the text file of the padded data.\n");
exit(EXIT_FAILURE);
}
printf("Task done - padded data can now be found in prj4input.txt \n\n");
outputptr=(double*)malloc(two_numdata*sizeof(double)); /*allocate space to the output complex array*/
if(outputptr==NULL)
{
printf("Error allocating memory to output array. The program will now terminate.\n");
exit(EXIT_FAILURE);
}
/*Call padme again for the output array*/
printf("The size of the output array will now be altered if necessary\n");
outputptr=padme(outputptr, &two_numdata);
numcomplex=two_numdata/2; /*Numcomplex is the number of complex numbers in the array, or half of
the array size in this case*/
/*Nested loops to determine whether a forward or inverse transform is to be carried out.
The loops use Y/N choices to allow the user to decide what they wish to do*/
while(isign==2)
{
printf("Do you wish to carry out a Fourier transform, y/n? ");
confirm=toupper(getchar()); /*Converts the input to upper case to catch, for example, both y and Y*/
if(confirm=='Y')
isign=1;
else if(confirm=='N')
{
/*If a forward transform is not required, check to see whether an inverse transform is*/
while(isign==2)
{
printf("Do you wish to carry out an inverse fourier transform, y/n? ");
confirm=toupper(getchar());
if(confirm=='Y')
isign=-1;
else if(confirm =='N')
isign=0;
}
}
}
printf("\n"); /*Puts an empty line between previous output and following output.
Makes final output easier to read and follow*/
if(isign==1)
{
printf("A Fourier transform of the data will now be carried out\n");
mycpxfft(spectrumptr, outputptr, numcomplex, isign);
}
else if(isign==-1)
{
printf("An inverse Fourier transform will now be carried out\n");
mycpxfft(spectrumptr, outputptr, numcomplex, isign);
}
else if(isign==0) /*If neither transform is required, the program will simply output the padded input array*/
printf("No transforming of the data will be carried out\n");
dummyptr=outputptr;
if ((sp=fopen(outfile,"w")) !=NULL)
{
for(i=1; i<=two_numdata; i++)
{
fprintf(sp,"%f \n",*dummyptr); /* NOTE: \n only puts in a line feed (chr = 10) not the
carriage-return linefeed pair (13 10) required by windows.
Excel/origin will still read*/
dummyptr++;
}
fclose(sp);
}
else
{
printf("Cannot write the text output file. \n");
exit(EXIT_FAILURE);
}
printf("Task done - all OK \n");
free(spectrumptr);
free(outputptr);
return(0);
}

>Any idea why the printf statements in my Y/N choice section would be executing twice? Do I need to clear the buffer?
I looked at the code in red, and it looks like the answer is yes to clearing the buffer. getchar() reads only the Y or N, but there is a newline also in the buffer. So you could either call getchar() one more time, or make a loop in case the user types "Yes". So the code in red contains two getchar() calls. After each you could add:

Code:

while (getchar() != '\n');

to clear the input stream. A second option would be to use fgets() to read a whole line of input.