annotate src/ext/kissfft/test/compfft.py @ 196:da283326bcd3 tip master

Update plugin versions in RDF
author Chris Cannam <cannam@all-day-breakfast.com>
date Fri, 28 Feb 2020 09:43:02 +0000
parents 5ed6e970541b
children
rev   line source
c@174 1 #!/usr/bin/env python
c@174 2
c@174 3 # use FFTPACK as a baseline
c@174 4 import FFT
c@174 5 from Numeric import *
c@174 6 import math
c@174 7 import random
c@174 8 import sys
c@174 9 import struct
c@174 10 import fft
c@174 11
c@174 12 pi=math.pi
c@174 13 e=math.e
c@174 14 j=complex(0,1)
c@174 15 lims=(-32768,32767)
c@174 16
c@174 17 def randbuf(n,cpx=1):
c@174 18 res = array( [ random.uniform( lims[0],lims[1] ) for i in range(n) ] )
c@174 19 if cpx:
c@174 20 res = res + j*randbuf(n,0)
c@174 21 return res
c@174 22
c@174 23 def main():
c@174 24 from getopt import getopt
c@174 25 import popen2
c@174 26 opts,args = getopt( sys.argv[1:],'u:n:Rt:' )
c@174 27 opts=dict(opts)
c@174 28 exitcode=0
c@174 29
c@174 30 util = opts.get('-u','./kf_float')
c@174 31
c@174 32 try:
c@174 33 dims = [ int(d) for d in opts['-n'].split(',')]
c@174 34 cpx = opts.get('-R') is None
c@174 35 fmt=opts.get('-t','f')
c@174 36 except KeyError:
c@174 37 sys.stderr.write("""
c@174 38 usage: compfft.py
c@174 39 -n d1[,d2,d3...] : FFT dimension(s)
c@174 40 -u utilname : see sample_code/fftutil.c, default = ./kf_float
c@174 41 -R : real-optimized version\n""")
c@174 42 sys.exit(1)
c@174 43
c@174 44 x = fft.make_random( dims )
c@174 45
c@174 46 cmd = '%s -n %s ' % ( util, ','.join([ str(d) for d in dims]) )
c@174 47 if cpx:
c@174 48 xout = FFT.fftnd(x)
c@174 49 xout = reshape(xout,(size(xout),))
c@174 50 else:
c@174 51 cmd += '-R '
c@174 52 xout = FFT.real_fft(x)
c@174 53
c@174 54 proc = popen2.Popen3( cmd , bufsize=len(x) )
c@174 55
c@174 56 proc.tochild.write( dopack( x , fmt ,cpx ) )
c@174 57 proc.tochild.close()
c@174 58 xoutcomp = dounpack( proc.fromchild.read( ) , fmt ,1 )
c@174 59 #xoutcomp = reshape( xoutcomp , dims )
c@174 60
c@174 61 sig = xout * conjugate(xout)
c@174 62 sigpow = sum( sig )
c@174 63
c@174 64 diff = xout-xoutcomp
c@174 65 noisepow = sum( diff * conjugate(diff) )
c@174 66
c@174 67 snr = 10 * math.log10(abs( sigpow / noisepow ) )
c@174 68 if snr<100:
c@174 69 print xout
c@174 70 print xoutcomp
c@174 71 exitcode=1
c@174 72 print 'NFFT=%s,SNR = %f dB' % (str(dims),snr)
c@174 73 sys.exit(exitcode)
c@174 74
c@174 75 def dopack(x,fmt,cpx):
c@174 76 x = reshape( x, ( size(x),) )
c@174 77 if cpx:
c@174 78 s = ''.join( [ struct.pack('ff',c.real,c.imag) for c in x ] )
c@174 79 else:
c@174 80 s = ''.join( [ struct.pack('f',c) for c in x ] )
c@174 81 return s
c@174 82
c@174 83 def dounpack(x,fmt,cpx):
c@174 84 uf = fmt * ( len(x) / 4 )
c@174 85 s = struct.unpack(uf,x)
c@174 86 if cpx:
c@174 87 return array(s[::2]) + array( s[1::2] )*j
c@174 88 else:
c@174 89 return array(s )
c@174 90
c@174 91 if __name__ == "__main__":
c@174 92 main()