Chris@184: #!/usr/bin/env python
Chris@184: 
Chris@184: from Numeric import *
Chris@184: from FFT import *
Chris@184: 
Chris@184: def make_random(len):
Chris@184:     import random
Chris@184:     res=[]
Chris@184:     for i in range(int(len)):
Chris@184:         r=random.uniform(-1,1)
Chris@184:         i=random.uniform(-1,1)
Chris@184:         res.append( complex(r,i) )
Chris@184:     return res
Chris@184: 
Chris@184: def slowfilter(sig,h):
Chris@184:     translen = len(h)-1
Chris@184:     return convolve(sig,h)[translen:-translen]
Chris@184: 
Chris@184: def nextpow2(x):
Chris@184:     return 2 ** math.ceil(math.log(x)/math.log(2))
Chris@184: 
Chris@184: def fastfilter(sig,h,nfft=None):
Chris@184:     if nfft is None:
Chris@184:         nfft = int( nextpow2( 2*len(h) ) )
Chris@184:     H = fft( h , nfft )
Chris@184:     scraplen = len(h)-1
Chris@184:     keeplen = nfft-scraplen
Chris@184:     res=[]
Chris@184:     isdone = 0
Chris@184:     lastidx = nfft
Chris@184:     idx0 = 0
Chris@184:     while not isdone:
Chris@184:         idx1 = idx0 + nfft
Chris@184:         if idx1 >= len(sig):
Chris@184:             idx1 = len(sig)
Chris@184:             lastidx = idx1-idx0
Chris@184:             if lastidx <= scraplen:
Chris@184:                 break
Chris@184:             isdone = 1
Chris@184:         Fss = fft(sig[idx0:idx1],nfft)
Chris@184:         fm = Fss * H
Chris@184:         m = inverse_fft(fm)
Chris@184:         res.append( m[scraplen:lastidx] )
Chris@184:         idx0 += keeplen
Chris@184:     return concatenate( res )
Chris@184: 
Chris@184: def main():
Chris@184:     import sys
Chris@184:     from getopt import getopt
Chris@184:     opts,args = getopt(sys.argv[1:],'rn:l:')
Chris@184:     opts=dict(opts)
Chris@184: 
Chris@184:     siglen = int(opts.get('-l',1e4 ) )
Chris@184:     hlen =50 
Chris@184:  
Chris@184:     nfft = int(opts.get('-n',128) )
Chris@184:     usereal = opts.has_key('-r')
Chris@184: 
Chris@184:     print 'nfft=%d'%nfft
Chris@184:     # make a signal
Chris@184:     sig = make_random( siglen )
Chris@184:     # make an impulse response
Chris@184:     h = make_random( hlen )
Chris@184:     #h=[1]*2+[0]*3
Chris@184:     if usereal:
Chris@184:         sig=[c.real for c in sig]
Chris@184:         h=[c.real for c in h]
Chris@184: 
Chris@184:     # perform MAC filtering
Chris@184:     yslow = slowfilter(sig,h)
Chris@184:     #print '<YSLOW>',yslow,'</YSLOW>'
Chris@184:     #yfast = fastfilter(sig,h,nfft)
Chris@184:     yfast = utilfastfilter(sig,h,nfft,usereal)
Chris@184:     #print yfast
Chris@184:     print 'len(yslow)=%d'%len(yslow)
Chris@184:     print 'len(yfast)=%d'%len(yfast)
Chris@184:     diff = yslow-yfast
Chris@184:     snr = 10*log10( abs( vdot(yslow,yslow) / vdot(diff,diff) ) )
Chris@184:     print 'snr=%s' % snr
Chris@184:     if snr < 10.0:
Chris@184:         print 'h=',h
Chris@184:         print 'sig=',sig[:5],'...'
Chris@184:         print 'yslow=',yslow[:5],'...'
Chris@184:         print 'yfast=',yfast[:5],'...'
Chris@184: 
Chris@184: def utilfastfilter(sig,h,nfft,usereal):
Chris@184:     import compfft
Chris@184:     import os
Chris@184:     open( 'sig.dat','w').write( compfft.dopack(sig,'f',not usereal) )
Chris@184:     open( 'h.dat','w').write( compfft.dopack(h,'f',not usereal) )
Chris@184:     if usereal: 
Chris@184:         util = './fastconvr' 
Chris@184:     else:
Chris@184:         util = './fastconv'
Chris@184:     cmd = 'time %s -n %d -i sig.dat -h h.dat -o out.dat' % (util, nfft)
Chris@184:     print cmd
Chris@184:     ec  = os.system(cmd)
Chris@184:     print 'exited->',ec
Chris@184:     return compfft.dounpack(open('out.dat').read(),'f',not usereal)
Chris@184: 
Chris@184: if __name__ == "__main__":
Chris@184:     main()