c@409: #!/usr/bin/env python
c@409: 
c@409: # use FFTPACK as a baseline
c@409: import FFT
c@409: from Numeric import *
c@409: import math
c@409: import random
c@409: import sys
c@409: import struct
c@409: import fft
c@409: 
c@409: pi=math.pi
c@409: e=math.e
c@409: j=complex(0,1)
c@409: lims=(-32768,32767)
c@409: 
c@409: def randbuf(n,cpx=1):
c@409:     res = array( [ random.uniform( lims[0],lims[1] ) for i in range(n) ] )
c@409:     if cpx:
c@409:         res = res + j*randbuf(n,0)
c@409:     return res
c@409: 
c@409: def main():
c@409:     from getopt import getopt
c@409:     import popen2
c@409:     opts,args = getopt( sys.argv[1:],'u:n:Rt:' )
c@409:     opts=dict(opts)
c@409:     exitcode=0
c@409: 
c@409:     util = opts.get('-u','./kf_float')
c@409: 
c@409:     try:
c@409:         dims = [ int(d) for d in opts['-n'].split(',')]
c@409:         cpx = opts.get('-R') is None
c@409:         fmt=opts.get('-t','f')
c@409:     except KeyError:
c@409:         sys.stderr.write("""
c@409:         usage: compfft.py 
c@409:         -n d1[,d2,d3...]  : FFT dimension(s)
c@409:         -u utilname : see sample_code/fftutil.c, default = ./kf_float
c@409:         -R : real-optimized version\n""")
c@409:         sys.exit(1)
c@409: 
c@409:     x = fft.make_random( dims )
c@409: 
c@409:     cmd = '%s -n %s ' % ( util, ','.join([ str(d) for d in dims]) )
c@409:     if cpx:
c@409:         xout = FFT.fftnd(x)
c@409:         xout = reshape(xout,(size(xout),))
c@409:     else:
c@409:         cmd += '-R '
c@409:         xout = FFT.real_fft(x)
c@409: 
c@409:     proc = popen2.Popen3( cmd , bufsize=len(x) )
c@409: 
c@409:     proc.tochild.write( dopack( x , fmt ,cpx ) )
c@409:     proc.tochild.close()
c@409:     xoutcomp = dounpack( proc.fromchild.read( ) , fmt ,1 )
c@409:     #xoutcomp = reshape( xoutcomp , dims )
c@409: 
c@409:     sig = xout * conjugate(xout)
c@409:     sigpow = sum( sig )
c@409: 
c@409:     diff = xout-xoutcomp
c@409:     noisepow = sum( diff * conjugate(diff) )
c@409: 
c@409:     snr = 10 * math.log10(abs( sigpow / noisepow ) )
c@409:     if snr<100:
c@409:         print xout
c@409:         print xoutcomp
c@409:         exitcode=1
c@409:     print 'NFFT=%s,SNR = %f dB' % (str(dims),snr)
c@409:     sys.exit(exitcode)
c@409: 
c@409: def dopack(x,fmt,cpx):
c@409:     x = reshape( x, ( size(x),) )
c@409:     if cpx:
c@409:         s = ''.join( [ struct.pack('ff',c.real,c.imag) for c in x ] )
c@409:     else:
c@409:         s = ''.join( [ struct.pack('f',c) for c in x ] )
c@409:     return s 
c@409: 
c@409: def dounpack(x,fmt,cpx):
c@409:     uf = fmt * ( len(x) / 4 )
c@409:     s = struct.unpack(uf,x)
c@409:     if cpx:
c@409:         return array(s[::2]) + array( s[1::2] )*j
c@409:     else:    
c@409:         return array(s )
c@409: 
c@409: if __name__ == "__main__":
c@409:     main()