Chris@184: #!/usr/bin/env python
Chris@184: 
Chris@184: import math
Chris@184: import sys
Chris@184: import os
Chris@184: import random
Chris@184: import struct
Chris@184: import popen2
Chris@184: import getopt
Chris@184: import numpy
Chris@184: 
Chris@184: pi=math.pi
Chris@184: e=math.e
Chris@184: j=complex(0,1)
Chris@184: 
Chris@184: doreal=0
Chris@184: 
Chris@184: datatype = os.environ.get('DATATYPE','float')
Chris@184: 
Chris@184: util = '../tools/fft_' + datatype
Chris@184: minsnr=90
Chris@184: if datatype == 'double':
Chris@184:     fmt='d'
Chris@184: elif datatype=='int16_t':
Chris@184:     fmt='h'
Chris@184:     minsnr=10
Chris@184: elif datatype=='int32_t':
Chris@184:     fmt='i'
Chris@184: elif datatype=='simd':
Chris@184:     fmt='4f'
Chris@184:     sys.stderr.write('testkiss.py does not yet test simd')
Chris@184:     sys.exit(0)
Chris@184: elif datatype=='float':
Chris@184:     fmt='f'
Chris@184: else:
Chris@184:     sys.stderr.write('unrecognized datatype %s\n' % datatype)
Chris@184:     sys.exit(1)
Chris@184:  
Chris@184: 
Chris@184: def dopack(x,cpx=1):
Chris@184:     x = numpy.reshape( x, ( numpy.size(x),) )
Chris@184:     
Chris@184:     if cpx:
Chris@184:         s = ''.join( [ struct.pack(fmt*2,c.real,c.imag) for c in x ] )
Chris@184:     else:
Chris@184:         s = ''.join( [ struct.pack(fmt,c.real) for c in x ] )
Chris@184:     return s
Chris@184: 
Chris@184: def dounpack(x,cpx):
Chris@184:     uf = fmt * ( len(x) / struct.calcsize(fmt) )
Chris@184:     s = struct.unpack(uf,x)
Chris@184:     if cpx:
Chris@184:         return numpy.array(s[::2]) + numpy.array( s[1::2] )*j
Chris@184:     else:
Chris@184:         return numpy.array(s )
Chris@184: 
Chris@184: def make_random(dims=[1]):
Chris@184:     res = []
Chris@184:     for i in range(dims[0]):
Chris@184:         if len(dims)==1:
Chris@184:             r=random.uniform(-1,1)
Chris@184:             if doreal:
Chris@184:                 res.append( r )
Chris@184:             else:
Chris@184:                 i=random.uniform(-1,1)
Chris@184:                 res.append( complex(r,i) )
Chris@184:         else:
Chris@184:             res.append( make_random( dims[1:] ) )
Chris@184:     return numpy.array(res)
Chris@184: 
Chris@184: def flatten(x):
Chris@184:     ntotal = numpy.size(x)
Chris@184:     return numpy.reshape(x,(ntotal,))
Chris@184: 
Chris@184: def randmat( ndims ):
Chris@184:     dims=[]
Chris@184:     for i in range( ndims ):
Chris@184:         curdim = int( random.uniform(2,5) )
Chris@184:         if doreal and i==(ndims-1):
Chris@184:             curdim = int(curdim/2)*2 # force even last dimension if real
Chris@184:         dims.append( curdim )
Chris@184:     return make_random(dims )
Chris@184: 
Chris@184: def test_fft(ndims):
Chris@184:     x=randmat( ndims )
Chris@184: 
Chris@184: 
Chris@184:     if doreal:
Chris@184:         xver = numpy.fft.rfftn(x)
Chris@184:     else:
Chris@184:         xver = numpy.fft.fftn(x)
Chris@184:     
Chris@184:     open('/tmp/fftexp.dat','w').write(dopack( flatten(xver) , True ) )
Chris@184: 
Chris@184:     x2=dofft(x,doreal)
Chris@184:     err = xver - x2
Chris@184:     errf = flatten(err)
Chris@184:     xverf = flatten(xver)
Chris@184:     errpow = numpy.vdot(errf,errf)+1e-10
Chris@184:     sigpow = numpy.vdot(xverf,xverf)+1e-10
Chris@184:     snr = 10*math.log10(abs(sigpow/errpow) )
Chris@184:     print 'SNR (compared to NumPy) : %.1fdB' % float(snr)
Chris@184: 
Chris@184:     if snr<minsnr:
Chris@184:         print 'xver=',xver
Chris@184:         print 'x2=',x2
Chris@184:         print 'err',err
Chris@184:         sys.exit(1)
Chris@184:  
Chris@184: def dofft(x,isreal):
Chris@184:     dims=list( numpy.shape(x) )
Chris@184:     x = flatten(x)
Chris@184: 
Chris@184:     scale=1
Chris@184:     if datatype=='int16_t':
Chris@184:         x = 32767 * x
Chris@184:         scale = len(x) / 32767.0
Chris@184:     elif datatype=='int32_t':
Chris@184:         x = 2147483647.0 * x
Chris@184:         scale = len(x) / 2147483647.0
Chris@184: 
Chris@184:     cmd='%s -n ' % util
Chris@184:     cmd += ','.join([str(d) for d in dims])
Chris@184:     if doreal:
Chris@184:         cmd += ' -R '
Chris@184: 
Chris@184:     print cmd
Chris@184:     p = popen2.Popen3(cmd )
Chris@184: 
Chris@184:     open('/tmp/fftin.dat','w').write(dopack( x , isreal==False ) )
Chris@184: 
Chris@184:     p.tochild.write( dopack( x , isreal==False ) )
Chris@184:     p.tochild.close()
Chris@184: 
Chris@184:     res = dounpack( p.fromchild.read() , 1 )
Chris@184:     open('/tmp/fftout.dat','w').write(dopack( flatten(res) , True ) )
Chris@184:     if doreal:
Chris@184:         dims[-1] = int( dims[-1]/2 ) + 1
Chris@184: 
Chris@184:     res = scale * res
Chris@184: 
Chris@184:     p.wait()
Chris@184:     return numpy.reshape(res,dims)
Chris@184: 
Chris@184: def main():
Chris@184:     opts,args = getopt.getopt(sys.argv[1:],'r')
Chris@184:     opts=dict(opts)
Chris@184: 
Chris@184:     global doreal
Chris@184:     doreal = opts.has_key('-r')
Chris@184: 
Chris@184:     if doreal:
Chris@184:         print 'Testing multi-dimensional real FFTs'
Chris@184:     else:
Chris@184:         print 'Testing multi-dimensional FFTs'
Chris@184: 
Chris@184:     for dim in range(1,4):
Chris@184:         test_fft( dim )
Chris@184: 
Chris@184: if __name__ == "__main__":
Chris@184:     main()
Chris@184: