1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
|
#!/usr/bin/env python
import getopt, sys
from cvxopt import matrix
from cvxopt.solvers import options
from cvxopt.modeling import op, variable, max
from math import cos, log10, pi
import pygtk
pygtk.require('2.0')
import gtk
import matplotlib
matplotlib.use('GTKAgg') # or 'GTK'
from matplotlib.backends.backend_gtkagg import FigureCanvasGTKAgg as FigureCanvas
import pylab
def frange(a,b,N):
return [ a+k*float((b-a))/N for k in range(N) ]
def design_lowpass(N, d1, wc, ws, solver=None, Q=50):
h = variable(N+1)
d1 = 10**(d1/20.0) # convert from dB
n1 = int(round(N*Q*wc/pi));
w1 = matrix(frange(0,wc,n1))
G1 = matrix([cos(wi*j) for j in range(N+1) for wi in w1], (n1,N+1))
n2 = int(round(N*Q*(pi-ws)/pi));
w2 = matrix(frange(ws,pi,n2))
G2 = matrix([cos(wi*j) for j in range(N+1) for wi in w2], (n2,N+1))
options['show_progress'] = 0
options['LPX_K_MSGLEV'] = 0
options['MSK_IPAR_LOG']= 0
op(max(abs(G2*h)), [G1*h <= d1, G1*h >= 1.0/d1]).solve(solver=solver)
return (h.value, max(abs(G2*h.value)))
def make_plot(h, d2, co, sb, pr, N, output=None):
w = matrix(frange(0,pi,N*50));
C = w*matrix(range(N+1), (1,N+1), 'd');
for i in range(len(C)): C[i] = cos(C[i])
fig = pylab.figure()
ax = fig.add_subplot(111)
ylim = [round((20*log10(d2)-40)/10)*10,10]
ax.plot(list(w/pi),[20*log10(abs(x)) for x in C*h])
ax.plot(2*[co],[-pr,ylim[0]],'g--')
ax.plot(2*[co],[10,pr],'g--')
ax.plot(2*[sb],[10,20*log10(d2)],'g--')
ax.plot([sb, 1],2*[20*log10(d2)],'g--')
ax.plot([0, co],2*[-pr],'g--')
ax.plot([0, co],2*[pr],'g--')
pylab.setp(ax,ylim=ylim)
pylab.setp(ax,xlabel="Normalized frequency")
pylab.setp(ax,ylabel="Attenuation [dB]")
ax.grid()
def usage():
print("""
Usage:
filterdemo_cli --cutoff=CO --stopband=SB --ripple=RP --order=N [options]
Arguments:
CO: normalized cutoff frequency
SB: normalized stopband frequency, 0.1 <= co < sb-0.1 <= 0.5,
RP: maximum passband ripple in dB, 0.01 <= rp <= 3,
N : filterorder, 5 <= or <= 50.
Options:
--solver = SOLVER One of default, mosek, glpk
--output = FILENAME Output filename.
""")
sys.exit(2)
def main():
try:
opts, args = getopt.getopt(sys.argv[1:], "",
['cutoff=', 'stopband=', 'ripple=', 'order=',
'solver=', 'output='])
except getopt.GetoptError:
usage()
if opts==[]: usage()
co, sb, pr, N, output = [None]*5
solver = "default"
try:
for o, a in opts:
if o == "--cutoff": co = float(a);
if o == "--stopband": sb = float(a);
if o == "--ripple": pr = float(a);
if o == "--order": N = int(a);
if o == "--solver": solver = a;
if o == "--output": output = a;
except: usage()
if None in [co, sb, pr, N]: usage()
if not (0.1 <= co < sb-0.01+1E-8 <= 0.5):
print("invalid cutoff and stopband frequencies")
usage()
if not (0.01 <= pr <= 3):
print("invalid value of passband ripple")
usage()
if not (5 <= N <= 50):
print("invalid filterorder")
usage()
if not solver in ['default','mosek','glpk']:
print("invalid solver")
usage()
try:
[h, d2] = design_lowpass(N, pr, co*pi, sb*pi, solver)
except:
print("Please tighten filter specifications.")
sys.exit(2)
make_plot(h, d2, co, sb, pr, N, output);
if (output != None): savefig(output)
pylab.show()
if __name__ == "__main__":
main()
|
