fixfft.cpp

这是我本人编写的一个32位定点小数运算的函数库。对于没有浮点运算器的场合

// FixFFT.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <FFT.h>
#include <conio.h>
#include <time.h>

#define N 1024
#define H ((double)(1<<14))

long fix_fft32(long fr[], long fi[], int m, int inverse);
extern "C" short fix_fft(short fr[], short fi[], int m, int inverse);

double randp()
{
	double r=(double)(rand()-RAND_MAX/2)/(RAND_MAX/2);
	r=(r>0)?r/4+0.75:r/4-0.75;
	return r;
}

int main()
{
	double datar[N],datai[N];
	short fix_datar[N],fix_datai[N];
	long fix_datar32[N],fix_datai32[N];

	CComplex<double> temp[N];
	short fix_temp_r[N],fix_temp_i[N];
	long fix_temp_r32[N],fix_temp_i32[N];


	srand((unsigned)time(NULL));
	double sum=0,sum32=0,maxd32=0;

	for (int ii=0;ii<1000;ii++) {
		if (ii%10==9) printf(".");
		for (int i=0; i<N; i++){
			datar[i]=randp(); datai[i]=randp();
	//		data[i]=cos(i*2*3.1415926535/N);
			fix_datar[i]=(short)(datar[i]*H); fix_datar32[i]=(long)(datar[i]*H*H);
			fix_datai[i]=(short)(datai[i]*H); fix_datai32[i]=(long)(datai[i]*H*H);

			temp[i].real()=datar[i]; temp[i].imag()=datai[i];
			fix_temp_r[i]=fix_datar[i]; fix_temp_i[i]=fix_datai[i];
			fix_temp_r32[i]=fix_datar32[i]; fix_temp_i32[i]=fix_datai32[i];
		}

		FFT(temp,N);
		for (i=0;i<N;i++) {
			temp[i]=temp[i]/(double)N;
			if (fabs(temp[i].real()*H*H)>INT_MAX || fabs(temp[i].imag()*H*H)>INT_MAX)
				printf("*");
		}
		//printf("\n");

		//iFFT(temp,N);

		fix_fft(fix_temp_r, fix_temp_i, 10, 0);
		fix_fft32(fix_temp_r32, fix_temp_i32, 10, 0);
		//fix_fft(fix_temp_r, fix_temp_i, 10, 1);
		//fix_fft32(fix_temp_r32, fix_temp_i32, 10, 1);

		for (i=0;i<N;i++) {
			double d=fabs(temp[i].real()-(double)fix_temp_r[i]/H)+fabs(temp[i].imag()-(double)fix_temp_i[i]/H);
			double d32=fabs(temp[i].real()-(double)fix_temp_r32[i]/H/H)+fabs(temp[i].imag()-(double)fix_temp_i32[i]/H/H);
			sum+=d; sum32+=d32; if (maxd32<d32) maxd32=d32;
		}
	}

	printf("\n");
	printf("delta=%lf(%lf)\tdelta32=%20.18lf(%lf)\n",sum/N/ii,sum/ii/N*H,sum32/ii/N,sum32/ii/N*H*H);
	printf("Max delta32=%20.18lf(%lf)\n",maxd32,maxd32*H*H);

	//for (i=0;i<N;i++) temp[i]=temp[i]*(double)N;
	//iFFT(temp,N);
	printf("The End.\n"); getch();
	return 0;
}