shun_iwasawa a35b8f
#!/usr/bin/env python
shun_iwasawa a35b8f
shun_iwasawa a35b8f
import FFT
shun_iwasawa a35b8f
import sys
shun_iwasawa a35b8f
import random
shun_iwasawa a35b8f
import re
shun_iwasawa a35b8f
j=complex(0,1)
shun_iwasawa a35b8f
shun_iwasawa a35b8f
def randvec(n,iscomplex):
shun_iwasawa a35b8f
    if iscomplex:
shun_iwasawa a35b8f
        return [
shun_iwasawa a35b8f
                int(random.uniform(-32768,32767) ) + j*int(random.uniform(-32768,32767) )
shun_iwasawa a35b8f
                for i in range(n) ]
shun_iwasawa a35b8f
    else:                
shun_iwasawa a35b8f
        return [ int(random.uniform(-32768,32767) ) for i in range(n) ]
shun_iwasawa a35b8f
    
shun_iwasawa a35b8f
def c_format(v,round=0):
shun_iwasawa a35b8f
    if round:
shun_iwasawa a35b8f
        return ','.join( [ '{%d,%d}' %(int(c.real),int(c.imag) ) for c in v ] ) 
shun_iwasawa a35b8f
    else:
shun_iwasawa a35b8f
        s= ','.join( [ '{%.60f ,%.60f }' %(c.real,c.imag) for c in v ] ) 
shun_iwasawa a35b8f
        return re.sub(r'\.?0+ ',' ',s)
shun_iwasawa a35b8f
shun_iwasawa a35b8f
def test_cpx( n,inverse ,short):
shun_iwasawa a35b8f
    v = randvec(n,1)
shun_iwasawa a35b8f
    scale = 1
shun_iwasawa a35b8f
    if short:
shun_iwasawa a35b8f
        minsnr=30
shun_iwasawa a35b8f
    else:
shun_iwasawa a35b8f
        minsnr=100
shun_iwasawa a35b8f
shun_iwasawa a35b8f
    if inverse:
shun_iwasawa a35b8f
        tvecout = FFT.inverse_fft(v)
shun_iwasawa a35b8f
        if short:
shun_iwasawa a35b8f
            scale = 1
shun_iwasawa a35b8f
        else:            
shun_iwasawa a35b8f
            scale = len(v)
shun_iwasawa a35b8f
    else:
shun_iwasawa a35b8f
        tvecout = FFT.fft(v)
shun_iwasawa a35b8f
        if short:
shun_iwasawa a35b8f
            scale = 1.0/len(v)
shun_iwasawa a35b8f
shun_iwasawa a35b8f
    tvecout = [ c * scale for c in tvecout ]
shun_iwasawa a35b8f
shun_iwasawa a35b8f
shun_iwasawa a35b8f
    s="""#define NFFT %d""" % len(v) + """
shun_iwasawa a35b8f
    {
shun_iwasawa a35b8f
        double snr;
shun_iwasawa a35b8f
        kiss_fft_cpx test_vec_in[NFFT] = { """  + c_format(v) + """};
shun_iwasawa a35b8f
        kiss_fft_cpx test_vec_out[NFFT] = {"""  + c_format( tvecout ) + """};
shun_iwasawa a35b8f
        kiss_fft_cpx testbuf[NFFT];
shun_iwasawa a35b8f
        void * cfg = kiss_fft_alloc(NFFT,%d,0,0);""" % inverse + """
shun_iwasawa a35b8f
shun_iwasawa a35b8f
        kiss_fft(cfg,test_vec_in,testbuf);
shun_iwasawa a35b8f
        snr = snr_compare(test_vec_out,testbuf,NFFT);
shun_iwasawa a35b8f
        printf("DATATYPE=" xstr(kiss_fft_scalar) ", FFT n=%d, inverse=%d, snr = %g dB\\n",NFFT,""" + str(inverse) + """,snr);
shun_iwasawa a35b8f
        if (snr<""" + str(minsnr) + """)
shun_iwasawa a35b8f
            exit_code++;
shun_iwasawa a35b8f
        free(cfg);
shun_iwasawa a35b8f
    }
shun_iwasawa a35b8f
#undef NFFT    
shun_iwasawa a35b8f
"""
shun_iwasawa a35b8f
    return s
shun_iwasawa a35b8f
shun_iwasawa a35b8f
def compare_func():
shun_iwasawa a35b8f
    s="""
shun_iwasawa a35b8f
#define xstr(s) str(s)
shun_iwasawa a35b8f
#define str(s) #s
shun_iwasawa a35b8f
double snr_compare( kiss_fft_cpx * test_vec_out,kiss_fft_cpx * testbuf, int n)
shun_iwasawa a35b8f
{
shun_iwasawa a35b8f
    int k;
shun_iwasawa a35b8f
    double sigpow,noisepow,err,snr,scale=0;
shun_iwasawa a35b8f
    kiss_fft_cpx err;
shun_iwasawa a35b8f
    sigpow = noisepow = .000000000000000000000000000001; 
shun_iwasawa a35b8f
shun_iwasawa a35b8f
    for (k=0;k
shun_iwasawa a35b8f
        sigpow += test_vec_out[k].r * test_vec_out[k].r + 
shun_iwasawa a35b8f
                  test_vec_out[k].i * test_vec_out[k].i;
shun_iwasawa a35b8f
        C_SUB(err,test_vec_out[k],testbuf[k].r);
shun_iwasawa a35b8f
        noisepow += err.r * err.r + err.i + err.i;
shun_iwasawa a35b8f
shun_iwasawa a35b8f
        if (test_vec_out[k].r)
shun_iwasawa a35b8f
            scale += testbuf[k].r / test_vec_out[k].r;
shun_iwasawa a35b8f
    }
shun_iwasawa a35b8f
    snr = 10*log10( sigpow / noisepow );
shun_iwasawa a35b8f
    scale /= n;
shun_iwasawa a35b8f
    if (snr<10)
shun_iwasawa a35b8f
        printf( "\\npoor snr, try a scaling factor %f\\n" , scale );
shun_iwasawa a35b8f
    return snr;
shun_iwasawa a35b8f
}
shun_iwasawa a35b8f
"""
shun_iwasawa a35b8f
    return s
shun_iwasawa a35b8f
shun_iwasawa a35b8f
def main():
shun_iwasawa a35b8f
shun_iwasawa a35b8f
    from getopt import getopt
shun_iwasawa a35b8f
    opts,args = getopt(sys.argv[1:],'s')
shun_iwasawa a35b8f
    opts = dict(opts)
shun_iwasawa a35b8f
    short = int( opts.has_key('-s') )
shun_iwasawa a35b8f
shun_iwasawa a35b8f
    fftsizes = args
shun_iwasawa a35b8f
    if not fftsizes:
shun_iwasawa a35b8f
        fftsizes = [ 1800 ]
shun_iwasawa a35b8f
    print '#include "kiss_fft.h"'
shun_iwasawa a35b8f
    print compare_func()
shun_iwasawa a35b8f
    print "int main() { int exit_code=0;\n"
shun_iwasawa a35b8f
    for n in fftsizes:
shun_iwasawa a35b8f
        n = int(n)
shun_iwasawa a35b8f
        print test_cpx(n,0,short)
shun_iwasawa a35b8f
        print test_cpx(n,1,short)
shun_iwasawa a35b8f
    print """
shun_iwasawa a35b8f
    return exit_code;
shun_iwasawa a35b8f
}
shun_iwasawa a35b8f
"""
shun_iwasawa a35b8f
shun_iwasawa a35b8f
if __name__ == "__main__":
shun_iwasawa a35b8f
    main()