Chris@366: #include <stdlib.h>
Chris@366: #include <string.h>
Chris@366: #include <stdio.h>
Chris@366: #include "kiss_fft.h"
Chris@366: #include "kiss_fftr.h"
Chris@366: #include <limits.h>
Chris@366: 
Chris@366: 
Chris@366: static
Chris@366: double two_tone_test( int nfft, int bin1,int bin2)
Chris@366: {
Chris@366:     kiss_fftr_cfg cfg	= NULL;
Chris@366:     kiss_fft_cpx *kout	= NULL;
Chris@366:     kiss_fft_scalar *tbuf = NULL;
Chris@366: 
Chris@366:     int i;
Chris@366:     double f1 = bin1*2*M_PI/nfft;
Chris@366:     double f2 = bin2*2*M_PI/nfft;
Chris@366:     double sigpow=0;
Chris@366:     double noisepow=0;
Chris@366: #if FIXED_POINT==32
Chris@366:     long maxrange = LONG_MAX;
Chris@366: #else
Chris@366:     long maxrange = SHRT_MAX;/* works fine for float too*/
Chris@366: #endif
Chris@366: 
Chris@366:     cfg = kiss_fftr_alloc(nfft , 0, NULL, NULL);
Chris@366:     tbuf	= KISS_FFT_MALLOC(nfft * sizeof(kiss_fft_scalar));
Chris@366:     kout	= KISS_FFT_MALLOC(nfft * sizeof(kiss_fft_cpx));
Chris@366: 
Chris@366:     /* generate a signal with two tones*/
Chris@366:     for (i = 0; i < nfft; i++) {
Chris@366: #ifdef USE_SIMD
Chris@366:         tbuf[i] = _mm_set1_ps( (maxrange>>1)*cos(f1*i)
Chris@366:                              + (maxrange>>1)*cos(f2*i) );
Chris@366: #else        
Chris@366:         tbuf[i] =  (maxrange>>1)*cos(f1*i)
Chris@366:                  + (maxrange>>1)*cos(f2*i);
Chris@366: #endif        
Chris@366:     }
Chris@366: 
Chris@366:     kiss_fftr(cfg, tbuf, kout);
Chris@366: 
Chris@366:     for (i=0;i < (nfft/2+1);++i) {
Chris@366: #ifdef USE_SIMD        
Chris@366:         double tmpr = (double)*(float*)&kout[i].r / (double)maxrange;
Chris@366:         double tmpi = (double)*(float*)&kout[i].i / (double)maxrange;
Chris@366: #else
Chris@366:         double tmpr = (double)kout[i].r / (double)maxrange;
Chris@366:         double tmpi = (double)kout[i].i / (double)maxrange;
Chris@366: #endif
Chris@366:         double mag2 = tmpr*tmpr + tmpi*tmpi;
Chris@366:         if (i!=0 && i!= nfft/2)
Chris@366:             mag2 *= 2; /* all bins except DC and Nyquist have symmetric counterparts implied*/
Chris@366: 
Chris@366:         /* if there is power in one of the expected bins, it is signal, otherwise noise*/
Chris@366:         if ( i!=bin1 && i != bin2 ) 
Chris@366:             noisepow += mag2;
Chris@366:         else
Chris@366:             sigpow += mag2;
Chris@366:     }
Chris@366:     kiss_fft_cleanup();
Chris@366:     /*printf("TEST %d,%d,%d noise @ %fdB\n",nfft,bin1,bin2,10*log10(noisepow/sigpow +1e-30) );*/
Chris@366:     return 10*log10(sigpow/(noisepow+1e-50) );
Chris@366: }
Chris@366: 
Chris@366: int main(int argc,char ** argv)
Chris@366: {
Chris@366:     int nfft = 4*2*2*3*5;
Chris@366:     if (argc>1) nfft = atoi(argv[1]);
Chris@366: 
Chris@366:     int i,j;
Chris@366:     double minsnr = 500;
Chris@366:     double maxsnr = -500;
Chris@366:     double snr;
Chris@366:     for (i=0;i<nfft/2;i+= (nfft>>4)+1) {
Chris@366:         for (j=i;j<nfft/2;j+=(nfft>>4)+7) {
Chris@366:             snr = two_tone_test(nfft,i,j);
Chris@366:             if (snr<minsnr) {
Chris@366:                 minsnr=snr;
Chris@366:             }
Chris@366:             if (snr>maxsnr) {
Chris@366:                 maxsnr=snr;
Chris@366:             }
Chris@366:         }
Chris@366:     }
Chris@366:     snr = two_tone_test(nfft,nfft/2,nfft/2);
Chris@366:     if (snr<minsnr) minsnr=snr;
Chris@366:     if (snr>maxsnr) maxsnr=snr;
Chris@366: 
Chris@366:     printf("TwoToneTest: snr ranges from %ddB to %ddB\n",(int)minsnr,(int)maxsnr);
Chris@366:     printf("sizeof(kiss_fft_scalar) = %d\n",(int)sizeof(kiss_fft_scalar) );
Chris@366:     return 0;
Chris@366: }