Chris@366: #include "kiss_fftr.h"
Chris@366: #include "_kiss_fft_guts.h"
Chris@366: #include <sys/times.h>
Chris@366: #include <time.h>
Chris@366: #include <unistd.h>
Chris@366: 
Chris@366: static double cputime(void)
Chris@366: {
Chris@366:     struct tms t;
Chris@366:     times(&t);
Chris@366:     return (double)(t.tms_utime + t.tms_stime)/  sysconf(_SC_CLK_TCK) ;
Chris@366: }
Chris@366: 
Chris@366: static
Chris@366: kiss_fft_scalar rand_scalar(void) 
Chris@366: {
Chris@366: #ifdef USE_SIMD
Chris@366:     return _mm_set1_ps(rand()-RAND_MAX/2);
Chris@366: #else
Chris@366:     kiss_fft_scalar s = (kiss_fft_scalar)(rand() -RAND_MAX/2);
Chris@366:     return s/2;
Chris@366: #endif
Chris@366: }
Chris@366: 
Chris@366: static
Chris@366: double snr_compare( kiss_fft_cpx * vec1,kiss_fft_cpx * vec2, int n)
Chris@366: {
Chris@366:     int k;
Chris@366:     double sigpow=1e-10,noisepow=1e-10,err,snr,scale=0;
Chris@366: 
Chris@366: #ifdef USE_SIMD
Chris@366:     float *fv1 = (float*)vec1;
Chris@366:     float *fv2 = (float*)vec2;
Chris@366:     for (k=0;k<8*n;++k) {
Chris@366:         sigpow += *fv1 * *fv1;
Chris@366:         err = *fv1 - *fv2;
Chris@366:         noisepow += err*err;
Chris@366:         ++fv1;
Chris@366:         ++fv2;
Chris@366:     }
Chris@366: #else
Chris@366:     for (k=0;k<n;++k) {
Chris@366:         sigpow += (double)vec1[k].r * (double)vec1[k].r + 
Chris@366:                   (double)vec1[k].i * (double)vec1[k].i;
Chris@366:         err = (double)vec1[k].r - (double)vec2[k].r;
Chris@366:         noisepow += err * err;
Chris@366:         err = (double)vec1[k].i - (double)vec2[k].i;
Chris@366:         noisepow += err * err;
Chris@366: 
Chris@366:         if (vec1[k].r)
Chris@366:             scale +=(double) vec2[k].r / (double)vec1[k].r;
Chris@366:     }
Chris@366: #endif    
Chris@366:     snr = 10*log10( sigpow / noisepow );
Chris@366:     scale /= n;
Chris@366:     if (snr<10) {
Chris@366:         printf( "\npoor snr, try a scaling factor %f\n" , scale );
Chris@366:         exit(1);
Chris@366:     }
Chris@366:     return snr;
Chris@366: }
Chris@366: 
Chris@366: #ifndef NUMFFTS
Chris@366: #define NUMFFTS 10000
Chris@366: #endif
Chris@366: 
Chris@366: 
Chris@366: int main(int argc,char ** argv)
Chris@366: {
Chris@366:     int nfft = 8*3*5;
Chris@366:     double ts,tfft,trfft;
Chris@366:     int i;
Chris@366:     if (argc>1)
Chris@366:         nfft = atoi(argv[1]);
Chris@366:     kiss_fft_cpx cin[nfft];
Chris@366:     kiss_fft_cpx cout[nfft];
Chris@366:     kiss_fft_cpx sout[nfft];
Chris@366:     kiss_fft_cfg  kiss_fft_state;
Chris@366:     kiss_fftr_cfg  kiss_fftr_state;
Chris@366: 
Chris@366:     kiss_fft_scalar rin[nfft+2];
Chris@366:     kiss_fft_scalar rout[nfft+2];
Chris@366:     kiss_fft_scalar zero;
Chris@366:     memset(&zero,0,sizeof(zero) ); // ugly way of setting short,int,float,double, or __m128 to zero
Chris@366: 
Chris@366:     srand(time(0));
Chris@366: 
Chris@366:     for (i=0;i<nfft;++i) {
Chris@366:         rin[i] = rand_scalar();
Chris@366:         cin[i].r = rin[i];
Chris@366:         cin[i].i = zero;
Chris@366:     }
Chris@366: 
Chris@366:     kiss_fft_state = kiss_fft_alloc(nfft,0,0,0);
Chris@366:     kiss_fftr_state = kiss_fftr_alloc(nfft,0,0,0);
Chris@366:     kiss_fft(kiss_fft_state,cin,cout);
Chris@366:     kiss_fftr(kiss_fftr_state,rin,sout);
Chris@366:     /*
Chris@366:     printf(" results from kiss_fft : (%f,%f), (%f,%f), (%f,%f) ...\n "
Chris@366:             , (float)cout[0].r , (float)cout[0].i
Chris@366:             , (float)cout[1].r , (float)cout[1].i
Chris@366:             , (float)cout[2].r , (float)cout[2].i); 
Chris@366:     printf(" results from kiss_fftr: (%f,%f), (%f,%f), (%f,%f) ...\n "
Chris@366:             , (float)sout[0].r , (float)sout[0].i
Chris@366:             , (float)sout[1].r , (float)sout[1].i
Chris@366:             , (float)sout[2].r , (float)sout[2].i); 
Chris@366:     */
Chris@366:         
Chris@366:     printf( "nfft=%d, inverse=%d, snr=%g\n",
Chris@366:             nfft,0, snr_compare(cout,sout,(nfft/2)+1) );
Chris@366:     ts = cputime();
Chris@366:     for (i=0;i<NUMFFTS;++i) {
Chris@366:         kiss_fft(kiss_fft_state,cin,cout);
Chris@366:     }
Chris@366:     tfft = cputime() - ts;
Chris@366:     
Chris@366:     ts = cputime();
Chris@366:     for (i=0;i<NUMFFTS;++i) {
Chris@366:         kiss_fftr( kiss_fftr_state, rin, cout );
Chris@366:         /* kiss_fftri(kiss_fftr_state,cout,rin); */
Chris@366:     }
Chris@366:     trfft = cputime() - ts;
Chris@366: 
Chris@366:     printf("%d complex ffts took %gs, real took %gs\n",NUMFFTS,tfft,trfft);
Chris@366: 
Chris@366:     free(kiss_fft_state);
Chris@366:     free(kiss_fftr_state);
Chris@366: 
Chris@366:     kiss_fft_state = kiss_fft_alloc(nfft,1,0,0);
Chris@366:     kiss_fftr_state = kiss_fftr_alloc(nfft,1,0,0);
Chris@366: 
Chris@366:     memset(cin,0,sizeof(cin));
Chris@366: #if 1
Chris@366:     for (i=1;i< nfft/2;++i) {
Chris@366:         //cin[i].r = (kiss_fft_scalar)(rand()-RAND_MAX/2);
Chris@366:         cin[i].r = rand_scalar();
Chris@366:         cin[i].i = rand_scalar();
Chris@366:     }
Chris@366: #else
Chris@366:     cin[0].r = 12000;
Chris@366:     cin[3].r = 12000;
Chris@366:     cin[nfft/2].r = 12000;
Chris@366: #endif
Chris@366: 
Chris@366:     // conjugate symmetry of real signal 
Chris@366:     for (i=1;i< nfft/2;++i) {
Chris@366:         cin[nfft-i].r = cin[i].r;
Chris@366:         cin[nfft-i].i = - cin[i].i;
Chris@366:     }
Chris@366: 
Chris@366:     kiss_fft(kiss_fft_state,cin,cout);
Chris@366:     kiss_fftri(kiss_fftr_state,cin,rout);
Chris@366:     /*
Chris@366:     printf(" results from inverse kiss_fft : (%f,%f), (%f,%f), (%f,%f), (%f,%f), (%f,%f) ...\n "
Chris@366:             , (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
Chris@366:             ); 
Chris@366: 
Chris@366:     printf(" results from inverse kiss_fftr: %f,%f,%f,%f,%f ... \n"
Chris@366:             ,(float)rout[0] ,(float)rout[1] ,(float)rout[2] ,(float)rout[3] ,(float)rout[4]);
Chris@366: */
Chris@366:     for (i=0;i<nfft;++i) {
Chris@366:         sout[i].r = rout[i];
Chris@366:         sout[i].i = zero;
Chris@366:     }
Chris@366: 
Chris@366:     printf( "nfft=%d, inverse=%d, snr=%g\n",
Chris@366:             nfft,1, snr_compare(cout,sout,nfft/2) );
Chris@366:     free(kiss_fft_state);
Chris@366:     free(kiss_fftr_state);
Chris@366: 
Chris@366:     return 0;
Chris@366: }