Thread: Working up to an FFT

  1. #1
    Registered User
    Join Date
    Oct 2007
    Posts
    16

    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:
    Code:
    /*prj04
      name: David Brown
      student number: 0500830
      usercode: phufai
    
    */
    
    #include<stdio.h>
    #include<stdlib.h>
    #include<math.h>
    #include <sys/types.h>
    #include <sys/stat.h>
    
    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 */
    
    	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;
    	
    	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*/
    		two_numdata=2*numdata;
    	
    		spectrumptr=(double*) malloc(two_numdata*sizeof(double));	/*request a pointer to enough
    										bytes for the complex array */
    		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<=numdata; i++) 
    			{			
    				fread(dummyptr,sizeof(float),1, fp);
    				printf("%g  |  ", *dummyptr);
    				dummyptr++;
    			}
    			
    			printf("\n");
    			fclose(fp);
    
    		}
    		else
    		{
    			printf("Cannot open data file. \n");
    			exit(0);
    		}
       	}
    	else
    	{
    		printf("Cannot read file attributes. \n");
    		exit(0);
    	}
    	
    	
    	
    	
    	
    	dummyptr=outputptr;
    
    	if ((sp=fopen(outfile,"w")) !=NULL)						/* write a text file as well */ 
    	{
    		for(i=1; i<=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);
    	return(0);
    }
    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?

  2. #2
    and the hat of int overfl Salem's Avatar
    Join Date
    Aug 2001
    Location
    The edge of the known universe
    Posts
    39,660
    > fread(dummyptr,sizeof(float),1, fp);
    fread won't perform the cast, so you need
    Code:
    float temp;
    fread(&temp,sizeof(float),1, fp);
    *dummyptr = temp;
    > dummyptr=outputptr;
    You want to reset this back to the start of the array, which is spectrumptr.
    outputptr is uninitialised junk.
    If you dance barefoot on the broken glass of undefined behaviour, you've got to expect the occasional cut.
    If at first you don't succeed, try writing your phone number on the exam paper.

  3. #3
    Registered User
    Join Date
    Oct 2007
    Posts
    16
    Thanks very much Salem, as usual you've helped me greatly.

    The reason for outputptr is that in the final version of the code, with the FFT added, it will be used. But at this stage it's not needed, which I hadn't twigged.

  4. #4
    Registered User
    Join Date
    Oct 2007
    Posts
    16

    Unhappy

    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.

    Any help would be greatly appreciated.
    lordbubonicus

  5. #5
    Registered User
    Join Date
    Oct 2007
    Posts
    16
    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.

    No progress on the second problem yet though.

  6. #6
    Registered User
    Join Date
    Oct 2001
    Posts
    2,934
    Here's some observations:
    Code:
    >double* padme(double *ptr, unsigned *long two_numdata)
    No need for the second parameter to be a pointer:
    Code:
    double* padme(double *ptr, unsigned long two_numdata)
    Code:
    >	*two_numdata=arraysize;
    You want to store two_numdata in arraysize, not the reverse:
    Code:
    	arraysize = two_numdata;
    Code:
    >	spectrumptr=padme(spectrumptr, &two_numdata);
    And the call would be:
    Code:
    	spectrumptr = padme(spectrumptr, two_numdata);

  7. #7
    Registered User
    Join Date
    Oct 2001
    Posts
    2,934
    And for you second problem:
    (1) isign is never initialized, so either initialize isign to some value other than 1, -1, or 0, or make your loop a do-while:
    Code:
    	do {
    	.
    	.
    	} while(isign!=1 && isign!=0 && isign!=-1);
    (2) confirm should be declared as an int, not a char, since getchar() returns an int.

    (3) toupper() is declared is <ctype.h>, so you should include this header.

  8. #8
    Registered User
    Join Date
    Oct 2007
    Posts
    16
    Thanks for the help Swoopy. However:

    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);
    }

  9. #9
    Registered User
    Join Date
    Oct 2007
    Posts
    16
    Ok, after some further error checking it seems as though the segmentation fault lies somewhere in this section of code:
    Code:
    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; 
    
    	}
    But I'm still stumped as to what's causing it.

  10. #10
    Registered User
    Join Date
    Oct 2007
    Posts
    16
    Update again, I've fixed the segmentation fault and my program appears to be working now. I say appears, becasue I'm currently still trying to test it.

  11. #11
    Registered User
    Join Date
    Oct 2001
    Posts
    2,934
    It looks like most of your for-loops go from 1 to two_numdata or 1-N, versus 0 to two_numdata-1 or 0 to N-1.

    >1: I'm still suffering from the warning about the signedness of the padme argument.
    off_t is likely defined as long (versus unsigned long) in sys/types.h

  12. #12
    Registered User
    Join Date
    Oct 2007
    Posts
    16
    Swoopy: The numbering of the for loops doesn't seem to have any effect if I change it. However, you're probably right that it should begin at 0 so I'll do so.

    Any idea why the printf statements in my Y/N choice section would be executing twice? Do I need to clear the buffer?

  13. #13
    Registered User
    Join Date
    Oct 2001
    Posts
    2,934
    >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.

  14. #14
    Registered User
    Join Date
    Oct 2007
    Posts
    16
    Thanks very much. I thought it might be something like that, but was unsure. You've really helped me a lot with this.

  15. #15
    Registered User
    Join Date
    Oct 2001
    Posts
    2,934
    Glad I could help and also glad to hear you fixed the segmentation fault.

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Function not working
    By sloopy in forum C Programming
    Replies: 31
    Last Post: 11-12-2005, 08:08 PM
  2. Program Not working Right
    By raven420smoke in forum C++ Programming
    Replies: 2
    Last Post: 09-16-2005, 03:21 AM
  3. Trying to eject D drive using code, but not working... :(
    By snowfrog in forum C++ Programming
    Replies: 3
    Last Post: 05-07-2005, 07:47 PM
  4. x on upper right corner not working
    By caduardo21 in forum Windows Programming
    Replies: 1
    Last Post: 02-20-2005, 08:35 PM
  5. cygwin -> unix , my code not working properly ;(
    By CyC|OpS in forum C Programming
    Replies: 4
    Last Post: 05-18-2002, 04:08 AM