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