52f71453d22bed0e2ad8ea1b21f2d399a6242833
[beacon-rx-433.git] / fir.py
1 # -*- coding: utf-8 -*-
2 """
3 Python code for basis FIR filter design
4 @author: Matti Pastell <matti.pastell@helsinki.fi>
5 http://mpastell.com
6 """
7
8 import sys
9 from pylab import *
10 import scipy.signal as signal
11
12 #Plot frequency and phase response
13 def mfreqz(b,a=1):
14     w,h = signal.freqz(b,a)
15     h_dB = 20 * log10 (abs(h))
16     subplot(211)
17     plot(w/max(w),h_dB)
18     ylim(-150, 5)
19     ylabel('Magnitude (db)')
20     xlabel(r'Normalized Frequency (x$\pi$rad/sample)')
21     title(r'Frequency response')
22     subplot(212)
23     h_Phase = unwrap(arctan2(imag(h),real(h)))
24     plot(w/max(w),h_Phase)
25     ylabel('Phase (radians)')
26     xlabel(r'Normalized Frequency (x$\pi$rad/sample)')
27     title(r'Phase response')
28     subplots_adjust(hspace=0.5)
29
30 #Plot step and impulse response
31 def impz(b,a=1):
32     l = len(b)
33     impulse = repeat(0.,l); impulse[0] =1.
34     x = arange(0,l)
35     response = signal.lfilter(b,a,impulse)
36     subplot(211)
37     stem(x, response)
38     ylabel('Amplitude')
39     xlabel(r'n (samples)')
40     title(r'Impulse response')
41     subplot(212)
42     step = cumsum(response)
43     stem(x, step)
44     ylabel('Amplitude')
45     xlabel(r'n (samples)')
46     title(r'Step response')
47     subplots_adjust(hspace=0.5)
48
49
50 Fe = 5000.
51 Fniq = Fe/2
52 transition = 1.001
53
54 n = 32
55 Ffond = 500
56 Fharm1 = 1500.
57 Fother = 850.
58
59 def bandpass(n, freq):
60     global Fe
61     a = signal.firwin(n, cutoff = (freq/Fniq)/transition, window = 'boxcar')
62     b = - signal.firwin(n, cutoff = (freq/Fniq)*transition, window = 'boxcar')
63     b[n/2] = b[n/2] + 1
64     d = - (a+b); d[n/2] = d[n/2] + 1
65     w,h = signal.freqz(d,1)
66     return d/max(abs(h))
67
68 def c_dump(fond, harm1, other):
69     f = fond * 32768
70     f = map(int, f)
71     h = harm1 * 32768
72     h = map(int, h)
73     o = other * 32768
74     o = map(int, o)
75     print "/* Apply the 3 filters on the bitfield. The results are stored in"
76     print " * fil_fond, fil_harm1 and fil_other. */"
77     print "static void apply_filters(void)"
78     print "{"
79     print "     fil_fond = 0;"
80     print "     fil_harm1 = 0;"
81     print "     fil_other = 0;"
82     print
83     for i in range(len(f)):
84         print " if (bitfield & (1UL << %d)) {"%(i)
85         print "         fil_fond += %d;"%(f[i])
86         print "         fil_harm1 += %d;"%(h[i])
87         print "         fil_other += %d;"%(o[i])
88         print " }"
89         print
90     print "#ifdef HOST_VERSION"
91     print "     fil_fond = saturate16(fil_fond);"
92     print "     fil_harm1 = saturate16(fil_harm1);"
93     print "     fil_other = saturate16(fil_other);"
94     print "#endif"
95     print "}"
96
97
98 figure(3)
99 fond = bandpass(n, Ffond)
100 harm1 = bandpass(n, Fharm1)
101 other = bandpass(n, Fother)
102 c_dump(fond, harm1, other)
103
104 #Frequency response
105 mfreqz(fond)
106 mfreqz(harm1)
107 mfreqz(other)
108
109 if len(sys.argv) >= 2:
110     savefig(sys.argv[1])
111 else:
112     show()