diff options
author | Jean-Marc Valin <Jean-Marc.Valin@csiro.au> | 2008-02-08 02:48:15 +0300 |
---|---|---|
committer | Jean-Marc Valin <Jean-Marc.Valin@csiro.au> | 2008-02-08 02:48:15 +0300 |
commit | e6586d21fad18d4cd4b72c8ff0f870db8ccbcbe6 (patch) | |
tree | 583d7c7bed9c4a000372bf28023f0e3bfa8a30c9 /tests/real-fft-test.c | |
parent | 4d0a7d0f1bc96afe8d1756f75c8b61645dd7e265 (diff) |
Real FFT cleanup, plus some testcases
Diffstat (limited to 'tests/real-fft-test.c')
-rw-r--r-- | tests/real-fft-test.c | 176 |
1 files changed, 176 insertions, 0 deletions
diff --git a/tests/real-fft-test.c b/tests/real-fft-test.c new file mode 100644 index 0000000..532ed8d --- /dev/null +++ b/tests/real-fft-test.c @@ -0,0 +1,176 @@ +#include "kiss_fftr.h" +#include "_kiss_fft_guts.h" +#include <sys/times.h> +#include <time.h> +#include <unistd.h> +#include <stdio.h> +#include <string.h> + +static double cputime(void) +{ + struct tms t; + times(&t); + return (double)(t.tms_utime + t.tms_stime)/ sysconf(_SC_CLK_TCK) ; +} + +static +kiss_fft_scalar rand_scalar(void) +{ +#ifdef USE_SIMD + return _mm_set1_ps(rand()-RAND_MAX/2); +#else + kiss_fft_scalar s = (kiss_fft_scalar)(rand() -RAND_MAX/2); + return s/2; +#endif +} + +static +double snr_compare( kiss_fft_cpx * vec1,kiss_fft_scalar * vec2, int n) +{ + int k; + double sigpow=1e-10, noisepow=1e-10, err,snr; + + for (k=1;k<n;++k) { + sigpow += (double)vec1[k].r * (double)vec1[k].r + + (double)vec1[k].i * (double)vec1[k].i; + err = (double)vec1[k].r - (double)vec2[2*k-1]; + noisepow += err * err; + err = (double)vec1[k].i - (double)vec2[2*k]; + noisepow += err * err; + + } + snr = 10*log10( sigpow / noisepow ); + if (snr<10) { + printf( "\npoor snr: %f\n", snr); + exit(1); + } + return snr; +} + +static +double snr_compare_scal( kiss_fft_scalar * vec1,kiss_fft_scalar * vec2, int n) +{ + int k; + double sigpow=1e-10, noisepow=1e-10, err,snr; + + for (k=1;k<n;++k) { + sigpow += (double)vec1[k] * (double)vec1[k]; + err = (double)vec1[k] - (double)vec2[k]; + noisepow += err * err; + } + snr = 10*log10( sigpow / noisepow ); + if (snr<10) { + printf( "\npoor snr: %f\n", snr); + exit(1); + } + return snr; +} +#define NFFT 8*3*5 + +#ifndef NUMFFTS +#define NUMFFTS 10000 +#endif + + +int main(void) +{ + double ts,tfft,trfft; + int i; + kiss_fft_cpx cin[NFFT]; + kiss_fft_cpx cout[NFFT]; + kiss_fft_scalar fin[NFFT]; + kiss_fft_scalar sout[NFFT]; + kiss_fft_cfg kiss_fft_state; + kiss_fftr_cfg kiss_fftr_state; + + kiss_fft_scalar rin[NFFT+2]; + kiss_fft_scalar rout[NFFT+2]; + kiss_fft_scalar zero; + memset(&zero,0,sizeof(zero) ); // ugly way of setting short,int,float,double, or __m128 to zero + + srand(time(0)); + + for (i=0;i<NFFT;++i) { + rin[i] = rand_scalar(); + cin[i].r = rin[i]; + cin[i].i = zero; + } + + kiss_fft_state = kiss_fft_alloc(NFFT,0,0,0); + kiss_fftr_state = kiss_fftr_alloc(NFFT,0,0,0); + kiss_fft(kiss_fft_state,cin,cout); + kiss_fftr(kiss_fftr_state,rin,sout); + + printf( "nfft=%d, inverse=%d, snr=%g\n", + NFFT,0, snr_compare(cout,sout,(NFFT/2)) ); + ts = cputime(); + for (i=0;i<NUMFFTS;++i) { + kiss_fft(kiss_fft_state,cin,cout); + } + tfft = cputime() - ts; + + ts = cputime(); + for (i=0;i<NUMFFTS;++i) { + kiss_fftr( kiss_fftr_state, rin, sout ); + /* kiss_fftri(kiss_fftr_state,cout,rin); */ + } + trfft = cputime() - ts; + + printf("%d complex ffts took %gs, real took %gs\n",NUMFFTS,tfft,trfft); + + free(kiss_fft_state); + free(kiss_fftr_state); + + kiss_fft_state = kiss_fft_alloc(NFFT,1,0,0); + kiss_fftr_state = kiss_fftr_alloc(NFFT,1,0,0); + + memset(cin,0,sizeof(cin)); +#if 1 + cin[0].r = rand_scalar(); + cin[NFFT/2].r = rand_scalar(); + for (i=1;i< NFFT/2;++i) { + //cin[i].r = (kiss_fft_scalar)(rand()-RAND_MAX/2); + cin[i].r = rand_scalar(); + cin[i].i = rand_scalar(); + } +#else + cin[0].r = 12000; + cin[3].r = 12000; + cin[NFFT/2].r = 12000; +#endif + + // conjugate symmetry of real signal + for (i=1;i< NFFT/2;++i) { + cin[NFFT-i].r = cin[i].r; + cin[NFFT-i].i = - cin[i].i; + } + + fin[0] = cin[0].r; + fin[NFFT-1] = cin[NFFT/2].r; + for (i=1;i< NFFT/2;++i) + { + fin[2*i-1] = cin[i].r; + fin[2*i] = cin[i].i; + } + + kiss_fft(kiss_fft_state,cin,cout); + kiss_fftri(kiss_fftr_state,fin,rout); + /* + printf(" results from inverse kiss_fft : (%f,%f), (%f,%f), (%f,%f), (%f,%f), (%f,%f) ...\n " + , (float)cout[0].r , (float)cout[0].i , (float)cout[1].r , (float)cout[1].i , (float)cout[2].r , (float)cout[2].i , (float)cout[3].r , (float)cout[3].i , (float)cout[4].r , (float)cout[4].i + ); + + printf(" results from inverse kiss_fftr: %f,%f,%f,%f,%f ... \n" + ,(float)rout[0] ,(float)rout[1] ,(float)rout[2] ,(float)rout[3] ,(float)rout[4]); +*/ + for (i=0;i<NFFT;++i) { + sout[i] = cout[i].r; + } + + printf( "nfft=%d, inverse=%d, snr=%g\n", + NFFT,1, snr_compare_scal(rout,sout,NFFT) ); + free(kiss_fft_state); + free(kiss_fftr_state); + + return 0; +} |