blob: 7e3f6bf1ab6154db5adaa6936a85454164698944 [file] [log] [blame]
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "fourier.h"
// RNG implemented localy to avoid library incongruences
#ifdef RAND_MAX
#undef RAND_MAX
#endif
#define RAND_MAX 32767
static unsigned long long int next = 1;
int rand( void ) {
next = next * 1103515245 + 12345;
return (unsigned int)(next / 65536) % RAND_MAX+1;
}
void srand( unsigned int seed ) {
next = seed;
}
// End of RNG implementation
/* Return 0 when V1 and V2 do not match within the allowed FP_ABSTOLERANCE. */
static inline int
check_FP(float V1, float V2) {
double AbsTolerance = FP_ABSTOLERANCE;
double Diff = fabs(V1 - V2);
if (Diff > AbsTolerance) {
fprintf(stderr, "A = %lf and B = %lf differ more than"
" FP_ABSTOLERANCE = %lf\n", V1, V2, AbsTolerance);
return 0;
}
return 1;
}
int main(int argc, char *argv[]) {
unsigned MAXSIZE;
unsigned MAXWAVES;
unsigned i,j;
float *RealIn;
float *ImagIn;
float *RealOut;
float *ImagOut;
float *RealOut_StrictFP;
float *ImagOut_StrictFP;
float *coeff;
float *amp;
int invfft=0;
if (argc<3)
{
printf("Usage: fft <waves> <length> -i\n");
printf("-i performs an inverse fft\n");
printf("make <waves> random sinusoids");
printf("<length> is the number of samples\n");
exit(-1);
}
else if (argc==4)
invfft = !strncmp(argv[3],"-i",2);
MAXSIZE=atoi(argv[2]);
MAXWAVES=atoi(argv[1]);
srand(1);
RealIn=(float*)malloc(sizeof(float)*MAXSIZE);
ImagIn=(float*)malloc(sizeof(float)*MAXSIZE);
RealOut=(float*)malloc(sizeof(float)*MAXSIZE);
ImagOut=(float*)malloc(sizeof(float)*MAXSIZE);
RealOut_StrictFP=(float*)malloc(sizeof(float)*MAXSIZE);
ImagOut_StrictFP=(float*)malloc(sizeof(float)*MAXSIZE);
coeff=(float*)malloc(sizeof(float)*MAXWAVES);
amp=(float*)malloc(sizeof(float)*MAXWAVES);
/* Makes MAXWAVES waves of random amplitude and period */
for(i=0;i<MAXWAVES;i++)
{
coeff[i] = i%1000;
amp[i] = i%1000;
}
for(i=0;i<MAXSIZE;i++)
{
/* RealIn[i]=rand();*/
RealIn[i]=0;
for(j=0;j<MAXWAVES;j++)
{
/* randomly select sin or cos */
if (rand()%2)
{
RealIn[i]+=coeff[j]*cos(amp[j]*i);
}
else
{
RealIn[i]+=coeff[j]*sin(amp[j]*i);
}
ImagIn[i]=0;
}
}
/* regular*/
fft_float (MAXSIZE,invfft,RealIn,ImagIn,RealOut,ImagOut);
fft_float_StrictFP (MAXSIZE,invfft,RealIn,ImagIn,RealOut_StrictFP,ImagOut_StrictFP);
printf("RealOut:\n");
for (i=0;i<MAXSIZE;i++) {
if (!check_FP(RealOut[i], RealOut_StrictFP[i]))
return 1;
printf("%f \t", RealOut_StrictFP[i]);
}
printf("\n");
printf("ImagOut:\n");
for (i=0;i<MAXSIZE;i++) {
if (!check_FP(ImagOut[i], ImagOut_StrictFP[i]))
return 1;
printf("%f \t", ImagOut_StrictFP[i]);
}
printf("\n");
free(RealIn);
free(ImagIn);
free(RealOut);
free(ImagOut);
free(RealOut_StrictFP);
free(ImagOut_StrictFP);
free(coeff);
free(amp);
exit(0);
}