fastfir.py
Upload User: ozl2332
Upload Date: 2009-12-28
Package Size: 38k
Code Size: 3k
Category:

Voice Compress

Development Platform:

C/C++

  1. #!/usr/bin/env python2.3
  2. from Numeric import *
  3. from FFT import *
  4. def make_random(len):
  5.     import random
  6.     res=[]
  7.     for i in range(int(len)):
  8.         r=random.uniform(-1,1)
  9.         i=random.uniform(-1,1)
  10.         res.append( complex(r,i) )
  11.     return res
  12. def slowfilter(sig,h):
  13.     translen = len(h)-1
  14.     return convolve(sig,h)[translen:-translen]
  15. def nextpow2(x):
  16.     return 2 ** math.ceil(math.log(x)/math.log(2))
  17. def fastfilter(sig,h,nfft=None):
  18.     if nfft is None:
  19.         nfft = int( nextpow2( 2*len(h) ) )
  20.     H = fft( h , nfft )
  21.     scraplen = len(h)-1
  22.     keeplen = nfft-scraplen
  23.     res=[]
  24.     isdone = 0
  25.     lastidx = nfft
  26.     idx0 = 0
  27.     while not isdone:
  28.         idx1 = idx0 + nfft
  29.         if idx1 >= len(sig):
  30.             idx1 = len(sig)
  31.             lastidx = idx1-idx0
  32.             if lastidx <= scraplen:
  33.                 break
  34.             isdone = 1
  35.         Fss = fft(sig[idx0:idx1],nfft)
  36.         fm = Fss * H
  37.         m = inverse_fft(fm)
  38.         res.append( m[scraplen:lastidx] )
  39.         idx0 += keeplen
  40.     return concatenate( res )
  41. def main():
  42.     import sys
  43.     from getopt import getopt
  44.     opts,args = getopt(sys.argv[1:],'rn:l:')
  45.     opts=dict(opts)
  46.     siglen = int(opts.get('-l',1e4 ) )
  47.     hlen =50 
  48.  
  49.     nfft = int(opts.get('-n',128) )
  50.     usereal = opts.has_key('-r')
  51.     print 'nfft=%d'%nfft
  52.     # make a signal
  53.     sig = make_random( siglen )
  54.     # make an impulse response
  55.     h = make_random( hlen )
  56.     #h=[1]*2+[0]*3
  57.     if usereal:
  58.         sig=[c.real for c in sig]
  59.         h=[c.real for c in h]
  60.     # perform MAC filtering
  61.     yslow = slowfilter(sig,h)
  62.     #print '<YSLOW>',yslow,'</YSLOW>'
  63.     #yfast = fastfilter(sig,h,nfft)
  64.     yfast = utilfastfilter(sig,h,nfft,usereal)
  65.     #print yfast
  66.     print 'len(yslow)=%d'%len(yslow)
  67.     print 'len(yfast)=%d'%len(yfast)
  68.     diff = yslow-yfast
  69.     snr = 10*log10( abs( vdot(yslow,yslow) / vdot(diff,diff) ) )
  70.     print 'snr=%s' % snr
  71.     if snr < 10.0:
  72.         print 'h=',h
  73.         print 'sig=',sig[:5],'...'
  74.         print 'yslow=',yslow[:5],'...'
  75.         print 'yfast=',yfast[:5],'...'
  76. def utilfastfilter(sig,h,nfft,usereal):
  77.     import compfft
  78.     import os
  79.     open( 'sig.dat','w').write( compfft.dopack(sig,'f',not usereal) )
  80.     open( 'h.dat','w').write( compfft.dopack(h,'f',not usereal) )
  81.     if usereal: 
  82.         util = './fastconvr' 
  83.     else:
  84.         util = './fastconv'
  85.     cmd = 'time %s -n %d -i sig.dat -h h.dat -o out.dat' % (util, nfft)
  86.     print cmd
  87.     ec  = os.system(cmd)
  88.     print 'exited->',ec
  89.     return compfft.dounpack(open('out.dat').read(),'f',not usereal)
  90. if __name__ == "__main__":
  91.     main()