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