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