annotate ext/kissfft/test/fastfir.py @ 409:1f1999b0f577

Bring in kissfft into this repo (formerly a subrepo, but the remote is not responding)

#!/usr/bin/env python

from Numeric import *
from FFT import *

def make_random(len):
    import random
    res=[]
    for i in range(int(len)):
        r=random.uniform(-1,1)
        i=random.uniform(-1,1)
        res.append( complex(r,i) )
    return res

def slowfilter(sig,h):
    translen = len(h)-1
    return convolve(sig,h)[translen:-translen]

def nextpow2(x):
    return 2 ** math.ceil(math.log(x)/math.log(2))

def fastfilter(sig,h,nfft=None):
    if nfft is None:
        nfft = int( nextpow2( 2*len(h) ) )
    H = fft( h , nfft )
    scraplen = len(h)-1
    keeplen = nfft-scraplen
    res=[]
    isdone = 0
    lastidx = nfft
    idx0 = 0
    while not isdone:
        idx1 = idx0 + nfft
        if idx1 >= len(sig):
            idx1 = len(sig)
            lastidx = idx1-idx0
            if lastidx <= scraplen:
                break
            isdone = 1
        Fss = fft(sig[idx0:idx1],nfft)
        fm = Fss * H
        m = inverse_fft(fm)
        res.append( m[scraplen:lastidx] )
        idx0 += keeplen
    return concatenate( res )

def main():
    import sys
    from getopt import getopt
    opts,args = getopt(sys.argv[1:],'rn:l:')
    opts=dict(opts)

    siglen = int(opts.get('-l',1e4 ) )
    hlen =50 
 
    nfft = int(opts.get('-n',128) )
    usereal = opts.has_key('-r')

    print 'nfft=%d'%nfft
    # make a signal
    sig = make_random( siglen )
    # make an impulse response
    h = make_random( hlen )
    #h=[1]*2+[0]*3
    if usereal:
        sig=[c.real for c in sig]
        h=[c.real for c in h]

    # perform MAC filtering
    yslow = slowfilter(sig,h)
    #print '<YSLOW>',yslow,'</YSLOW>'
    #yfast = fastfilter(sig,h,nfft)
    yfast = utilfastfilter(sig,h,nfft,usereal)
    #print yfast
    print 'len(yslow)=%d'%len(yslow)
    print 'len(yfast)=%d'%len(yfast)
    diff = yslow-yfast
    snr = 10*log10( abs( vdot(yslow,yslow) / vdot(diff,diff) ) )
    print 'snr=%s' % snr
    if snr < 10.0:
        print 'h=',h
        print 'sig=',sig[:5],'...'
        print 'yslow=',yslow[:5],'...'
        print 'yfast=',yfast[:5],'...'

def utilfastfilter(sig,h,nfft,usereal):
    import compfft
    import os
    open( 'sig.dat','w').write( compfft.dopack(sig,'f',not usereal) )
    open( 'h.dat','w').write( compfft.dopack(h,'f',not usereal) )
    if usereal: 
        util = './fastconvr' 
    else:
        util = './fastconv'
    cmd = 'time %s -n %d -i sig.dat -h h.dat -o out.dat' % (util, nfft)
    print cmd
    ec  = os.system(cmd)
    print 'exited->',ec
    return compfft.dounpack(open('out.dat').read(),'f',not usereal)

if __name__ == "__main__":
    main()