103 lines
2.6 KiB
Python
103 lines
2.6 KiB
Python
|
#!/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()
|