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