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