|
Packit |
664db3 |
#ifdef HAVE_CONFIG_H
|
|
Packit |
664db3 |
#include "config.h"
|
|
Packit |
664db3 |
#endif
|
|
Packit |
664db3 |
|
|
Packit |
664db3 |
#include "kiss_fftr.c"
|
|
Packit |
664db3 |
#include "kiss_fft.c"
|
|
Packit |
664db3 |
#include "_kiss_fft_guts.h"
|
|
Packit |
664db3 |
#include <stdio.h>
|
|
Packit |
664db3 |
#include <string.h>
|
|
Packit |
664db3 |
|
|
Packit |
664db3 |
int ret=0;
|
|
Packit |
664db3 |
|
|
Packit |
664db3 |
static
|
|
Packit |
664db3 |
kiss_fft_scalar rand_scalar(void)
|
|
Packit |
664db3 |
{
|
|
Packit |
664db3 |
return (rand()%32767)-16384;
|
|
Packit |
664db3 |
}
|
|
Packit |
664db3 |
|
|
Packit |
664db3 |
static
|
|
Packit |
664db3 |
double snr_compare( kiss_fft_cpx * vec1,kiss_fft_scalar * vec2, int n)
|
|
Packit |
664db3 |
{
|
|
Packit |
664db3 |
int k;
|
|
Packit |
664db3 |
double sigpow=1e-10, noisepow=1e-10, err,snr;
|
|
Packit |
664db3 |
|
|
Packit |
664db3 |
vec1[0].i = vec1[n].r;
|
|
Packit |
664db3 |
for (k=0;k
|
|
Packit |
664db3 |
sigpow += (double)vec1[k].r * (double)vec1[k].r +
|
|
Packit |
664db3 |
(double)vec1[k].i * (double)vec1[k].i;
|
|
Packit |
664db3 |
err = (double)vec1[k].r - (double)vec2[2*k];
|
|
Packit |
664db3 |
/*printf ("%f %f\n", (double)vec1[k].r, (double)vec2[2*k]);*/
|
|
Packit |
664db3 |
noisepow += err * err;
|
|
Packit |
664db3 |
err = (double)vec1[k].i - (double)vec2[2*k+1];
|
|
Packit |
664db3 |
/*printf ("%f %f\n", (double)vec1[k].i, (double)vec2[2*k+1]);*/
|
|
Packit |
664db3 |
noisepow += err * err;
|
|
Packit |
664db3 |
|
|
Packit |
664db3 |
}
|
|
Packit |
664db3 |
snr = 10*log10( sigpow / noisepow );
|
|
Packit |
664db3 |
if (snr<60) {
|
|
Packit |
664db3 |
printf( "** poor snr: %f **\n", snr);
|
|
Packit |
664db3 |
ret = 1;
|
|
Packit |
664db3 |
}
|
|
Packit |
664db3 |
return snr;
|
|
Packit |
664db3 |
}
|
|
Packit |
664db3 |
|
|
Packit |
664db3 |
static
|
|
Packit |
664db3 |
double snr_compare_scal( kiss_fft_scalar * vec1,kiss_fft_scalar * vec2, int n)
|
|
Packit |
664db3 |
{
|
|
Packit |
664db3 |
int k;
|
|
Packit |
664db3 |
double sigpow=1e-10, noisepow=1e-10, err,snr;
|
|
Packit |
664db3 |
|
|
Packit |
664db3 |
for (k=0;k
|
|
Packit |
664db3 |
sigpow += (double)vec1[k] * (double)vec1[k];
|
|
Packit |
664db3 |
err = (double)vec1[k] - (double)vec2[k];
|
|
Packit |
664db3 |
noisepow += err * err;
|
|
Packit |
664db3 |
}
|
|
Packit |
664db3 |
snr = 10*log10( sigpow / noisepow );
|
|
Packit |
664db3 |
if (snr<60) {
|
|
Packit |
664db3 |
printf( "\npoor snr: %f\n", snr);
|
|
Packit |
664db3 |
ret = 1;
|
|
Packit |
664db3 |
}
|
|
Packit |
664db3 |
return snr;
|
|
Packit |
664db3 |
}
|
|
Packit |
664db3 |
#ifdef RADIX_TWO_ONLY
|
|
Packit |
664db3 |
#define NFFT 1024
|
|
Packit |
664db3 |
#else
|
|
Packit |
664db3 |
#define NFFT 8*3*5
|
|
Packit |
664db3 |
#endif
|
|
Packit |
664db3 |
|
|
Packit |
664db3 |
#ifndef NUMFFTS
|
|
Packit |
664db3 |
#define NUMFFTS 10000
|
|
Packit |
664db3 |
#endif
|
|
Packit |
664db3 |
|
|
Packit |
664db3 |
|
|
Packit |
664db3 |
int main(void)
|
|
Packit |
664db3 |
{
|
|
Packit |
664db3 |
int i;
|
|
Packit |
664db3 |
kiss_fft_cpx cin[NFFT];
|
|
Packit |
664db3 |
kiss_fft_cpx cout[NFFT];
|
|
Packit |
664db3 |
kiss_fft_scalar fin[NFFT];
|
|
Packit |
664db3 |
kiss_fft_scalar sout[NFFT];
|
|
Packit |
664db3 |
kiss_fft_cfg kiss_fft_state;
|
|
Packit |
664db3 |
kiss_fftr_cfg kiss_fftr_state;
|
|
Packit |
664db3 |
|
|
Packit |
664db3 |
kiss_fft_scalar rin[NFFT+2];
|
|
Packit |
664db3 |
kiss_fft_scalar rout[NFFT+2];
|
|
Packit |
664db3 |
kiss_fft_scalar zero;
|
|
Packit |
664db3 |
memset(&zero,0,sizeof(zero) ); // ugly way of setting short,int,float,double, or __m128 to zero
|
|
Packit |
664db3 |
|
|
Packit |
664db3 |
for (i=0;i
|
|
Packit |
664db3 |
rin[i] = rand_scalar();
|
|
Packit |
664db3 |
#if defined(FIXED_POINT) && defined(DOUBLE_PRECISION)
|
|
Packit |
664db3 |
rin[i] *= 32768;
|
|
Packit |
664db3 |
#endif
|
|
Packit |
664db3 |
cin[i].r = rin[i];
|
|
Packit |
664db3 |
cin[i].i = zero;
|
|
Packit |
664db3 |
}
|
|
Packit |
664db3 |
|
|
Packit |
664db3 |
kiss_fft_state = kiss_fft_alloc(NFFT,0,0);
|
|
Packit |
664db3 |
kiss_fftr_state = kiss_fftr_alloc(NFFT,0,0);
|
|
Packit |
664db3 |
kiss_fft(kiss_fft_state,cin,cout);
|
|
Packit |
664db3 |
kiss_fftr(kiss_fftr_state,rin,sout);
|
|
Packit |
664db3 |
|
|
Packit |
664db3 |
printf( "nfft=%d, inverse=%d, snr=%g\n",
|
|
Packit |
664db3 |
NFFT,0, snr_compare(cout,sout,(NFFT/2)) );
|
|
Packit |
664db3 |
|
|
Packit |
664db3 |
memset(cin,0,sizeof(cin));
|
|
Packit |
664db3 |
cin[0].r = rand_scalar();
|
|
Packit |
664db3 |
cin[NFFT/2].r = rand_scalar();
|
|
Packit |
664db3 |
for (i=1;i< NFFT/2;++i) {
|
|
Packit |
664db3 |
//cin[i].r = (kiss_fft_scalar)(rand()-RAND_MAX/2);
|
|
Packit |
664db3 |
cin[i].r = rand_scalar();
|
|
Packit |
664db3 |
cin[i].i = rand_scalar();
|
|
Packit |
664db3 |
}
|
|
Packit |
664db3 |
|
|
Packit |
664db3 |
// conjugate symmetry of real signal
|
|
Packit |
664db3 |
for (i=1;i< NFFT/2;++i) {
|
|
Packit |
664db3 |
cin[NFFT-i].r = cin[i].r;
|
|
Packit |
664db3 |
cin[NFFT-i].i = - cin[i].i;
|
|
Packit |
664db3 |
}
|
|
Packit |
664db3 |
|
|
Packit |
664db3 |
|
|
Packit |
664db3 |
#ifdef FIXED_POINT
|
|
Packit |
664db3 |
#ifdef DOUBLE_PRECISION
|
|
Packit |
664db3 |
for (i=0;i< NFFT;++i) {
|
|
Packit |
664db3 |
cin[i].r *= 32768;
|
|
Packit |
664db3 |
cin[i].i *= 32768;
|
|
Packit |
664db3 |
}
|
|
Packit |
664db3 |
#endif
|
|
Packit |
664db3 |
for (i=0;i< NFFT;++i) {
|
|
Packit |
664db3 |
cin[i].r /= NFFT;
|
|
Packit |
664db3 |
cin[i].i /= NFFT;
|
|
Packit |
664db3 |
}
|
|
Packit |
664db3 |
#endif
|
|
Packit |
664db3 |
|
|
Packit |
664db3 |
fin[0] = cin[0].r;
|
|
Packit |
664db3 |
fin[1] = cin[NFFT/2].r;
|
|
Packit |
664db3 |
for (i=1;i< NFFT/2;++i)
|
|
Packit |
664db3 |
{
|
|
Packit |
664db3 |
fin[2*i] = cin[i].r;
|
|
Packit |
664db3 |
fin[2*i+1] = cin[i].i;
|
|
Packit |
664db3 |
}
|
|
Packit |
664db3 |
|
|
Packit |
664db3 |
kiss_ifft(kiss_fft_state,cin,cout);
|
|
Packit |
664db3 |
kiss_fftri(kiss_fftr_state,fin,rout);
|
|
Packit |
664db3 |
/*
|
|
Packit |
664db3 |
printf(" results from inverse kiss_fft : (%f,%f), (%f,%f), (%f,%f), (%f,%f), (%f,%f) ...\n "
|
|
Packit |
664db3 |
, (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
|
|
Packit |
664db3 |
);
|
|
Packit |
664db3 |
|
|
Packit |
664db3 |
printf(" results from inverse kiss_fftr: %f,%f,%f,%f,%f ... \n"
|
|
Packit |
664db3 |
,(float)rout[0] ,(float)rout[1] ,(float)rout[2] ,(float)rout[3] ,(float)rout[4]);
|
|
Packit |
664db3 |
*/
|
|
Packit |
664db3 |
for (i=0;i
|
|
Packit |
664db3 |
sout[i] = cout[i].r;
|
|
Packit |
664db3 |
}
|
|
Packit |
664db3 |
|
|
Packit |
664db3 |
printf( "nfft=%d, inverse=%d, snr=%g\n",
|
|
Packit |
664db3 |
NFFT,1, snr_compare_scal(rout,sout,NFFT) );
|
|
Packit |
664db3 |
free(kiss_fft_state);
|
|
Packit |
664db3 |
free(kiss_fftr_state);
|
|
Packit |
664db3 |
|
|
Packit |
664db3 |
return ret;
|
|
Packit |
664db3 |
}
|