← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to Lib/signal/filter_design.py

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
Import upstream version 0.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/signal/filter_design.py
 
1
"""Filter design.
 
2
"""
 
3
 
 
4
import scipy_base
 
5
from scipy_base.fastumath import *
 
6
from scipy_base import atleast_1d, poly, polyval, roots, imag, real, asarray,\
 
7
     allclose, Float, resize, pi, concatenate, absolute, logspace, c_
 
8
from scipy import comb, special, optimize, linalg
 
9
import string, types
 
10
 
 
11
 
 
12
MLab = scipy_base
 
13
Num = MLab
 
14
abs = absolute
 
15
 
 
16
def findfreqs(num, den, N):
 
17
    ep = atleast_1d(roots(den))+0j
 
18
    tz = atleast_1d(roots(num))+0j
 
19
 
 
20
    if len(ep) == 0:
 
21
        ep = atleast_1d(-1000)+0j
 
22
 
 
23
    ez = c_[Num.compress(ep.imag >=0, ep), Num.compress((abs(tz) < 1e5) & (tz.imag >=0),tz)]
 
24
 
 
25
    integ = abs(ez) < 1e-10
 
26
    hfreq = Num.around(Num.log10(Num.max(3*abs(ez.real + integ)+1.5*ez.imag))+0.5)
 
27
    lfreq = Num.around(Num.log10(0.1*Num.min(abs(real(ez+integ))+2*ez.imag))-0.5)
 
28
 
 
29
    w = logspace(lfreq, hfreq, N)
 
30
    return w
 
31
 
 
32
def freqs(b,a,worN=None,plot=None):
 
33
    """Compute frequency response of analog filter.
 
34
 
 
35
    Description:
 
36
 
 
37
       Given the numerator (b) and denominator (a) of a filter compute its
 
38
       frequency response.
 
39
 
 
40
               b[0]*(jw)**(nb-1) + b[1]*(jw)**(nb-2) + ... + b[nb-1]
 
41
       H(w) = --------------------------------------------------------
 
42
               a[0]*(jw)**(na-1) + a[1]*(jw)**(na-2) + ... + a[na-1]
 
43
 
 
44
    Inputs:
 
45
 
 
46
       b, a --- the numerator and denominator of a linear filter.
 
47
       worN --- If None, then compute at 200 frequencies around the interesting
 
48
                parts of the response curve (determined by pole-zero locations).
 
49
                If a single integer, the compute at that many frequencies.
 
50
                Otherwise, compute the response at frequencies given in worN.
 
51
    Outputs: (w,h)
 
52
 
 
53
       w -- The frequencies at which h was computed.
 
54
       h -- The frequency response.
 
55
    """
 
56
    if worN is None:
 
57
        w = findfreqs(b,a,200)
 
58
    elif isinstance(worN, types.IntType):
 
59
        N = worN
 
60
        w = findfreqs(b,a,N)
 
61
    else:
 
62
        w = worN
 
63
    w = atleast_1d(w)
 
64
    s = 1j*w
 
65
    h = polyval(b, s) / polyval(a, s)
 
66
    if not plot is None:
 
67
        plot(w, h)
 
68
    return w, h
 
69
 
 
70
def freqz(b, a, worN=None, whole=0, plot=None):
 
71
    """Compute frequency response of a digital filter.
 
72
 
 
73
    Description:
 
74
 
 
75
       Given the numerator (b) and denominator (a) of a digital filter compute
 
76
       its frequency response.
 
77
 
 
78
                  jw               -jw            -jmw
 
79
           jw  B(e)    b[0] + b[1]e + .... + b[m]e
 
80
        H(e) = ---- = ------------------------------------
 
81
                  jw               -jw            -jnw
 
82
               A(e)    a[0] + a[2]e + .... + a[n]e             
 
83
 
 
84
    Inputs:
 
85
 
 
86
       b, a --- the numerator and denominator of a linear filter.
 
87
       worN --- If None, then compute at 512 frequencies around the unit circle.
 
88
                If a single integer, the compute at that many frequencies.
 
89
                Otherwise, compute the response at frequencies given in worN
 
90
       whole -- Normally, frequencies are computed from 0 to pi (upper-half of
 
91
                unit-circle.  If whole is non-zero compute frequencies from 0
 
92
                to 2*pi.
 
93
 
 
94
    Outputs: (w,h)
 
95
 
 
96
       w -- The frequencies at which h was computed.
 
97
       h -- The frequency response.
 
98
    """
 
99
    b, a = map(atleast_1d, (b,a))
 
100
    if whole:
 
101
        lastpoint = 2*pi
 
102
    else:
 
103
        lastpoint = pi
 
104
    if worN is None:
 
105
        N = 512
 
106
        w = Num.arange(0,lastpoint,lastpoint/N)
 
107
    elif isinstance(worN, types.IntType):
 
108
        N = worN
 
109
        w = Num.arange(0,lastpoint,lastpoint/N)
 
110
    else:
 
111
        w = worN
 
112
    w = atleast_1d(w)
 
113
    zm1 = exp(-1j*w)
 
114
    h = polyval(b[::-1], zm1) / polyval(a[::-1], zm1)
 
115
    if not plot is None:
 
116
        plot(w, h)
 
117
    return w, h
 
118
    
 
119
def tf2zpk(b,a):
 
120
    """Return zero, pole, gain (z,p,k) representation from a numerator,
 
121
    denominator representation of a linear filter.
 
122
    """
 
123
    b,a = normalize(b,a)
 
124
    b = (b+0.0) / a[0]
 
125
    a = (a+0.0) / a[0]
 
126
    k = b[0]
 
127
    b /= b[0]
 
128
    z = roots(b)
 
129
    p = roots(a)
 
130
    return z, p, k
 
131
    
 
132
def zpk2tf(z,p,k):
 
133
    """Return polynomial transfer function representation from zeros
 
134
    and poles
 
135
 
 
136
    Inputs:
 
137
 
 
138
      z, p --- sequences representing the zeros and poles.
 
139
      k --- system gain.
 
140
 
 
141
    Outputs: (b,a)
 
142
 
 
143
      b, a --- numerator and denominator polynomials.
 
144
    """
 
145
    z = atleast_1d(z)
 
146
    k = atleast_1d(k)
 
147
    if len(z.shape) > 1:
 
148
        temp = poly(z[0])
 
149
        b = zeros((z.shape[0], z.shape[1]+1), temp.typecode())
 
150
        if len(k) == 1:
 
151
            k = [k[0]]*z.shape[0]
 
152
        for i in range(z.shape[0]):
 
153
            b[i] = k[i] * poly(z[i])
 
154
    else:
 
155
        b = k * poly(z)
 
156
    a = poly(p)
 
157
    return b, a
 
158
 
 
159
def normalize(b,a):
 
160
    """Normalize polynomial representation of a transfer function.
 
161
    """
 
162
    b,a = map(atleast_1d,(b,a))
 
163
    if len(a.shape) != 1:
 
164
        raise ValueError, "Denominator polynomial must be rank-1 array."
 
165
    if len(b.shape) > 2:
 
166
        raise ValueError, "Numerator polynomial must be rank-1 or rank-2 array."
 
167
    if len(b.shape) == 1:
 
168
        b = asarray([b],b.typecode())
 
169
    while a[0] == 0.0 and len(a) > 1:
 
170
        a = a[1:]
 
171
    while allclose(b[:,0], 0, rtol=1e-14) and (b.shape[-1] > 1):
 
172
        b = b[:,1:]
 
173
    if b.shape[0] == 1:
 
174
        b = b[0]        
 
175
    outb = b * (1.0) / a[0]
 
176
    outa = a * (1.0) / a[0]
 
177
    return outb, outa
 
178
 
 
179
 
 
180
def lp2lp(b,a,wo=1.0):
 
181
    """Return a low-pass filter with cuttoff frequency wo
 
182
    from a low-pass filter prototype with unity cutoff frequency.
 
183
    """
 
184
    a,b = map(atleast_1d,(a,b))
 
185
    if type(wo) is type(a):
 
186
        wo = wo[0]
 
187
    wo = float(wo)
 
188
    d = len(a)
 
189
    n = len(b)
 
190
    M = max((d,n))
 
191
    pwo = pow(wo,Num.arange(M-1,-1,-1))
 
192
    start1 = max((n-d,0))
 
193
    start2 = max((d-n,0))
 
194
    b = b * pwo[start1]/pwo[start2:]
 
195
    a = a * pwo[start1]/pwo[start1:]
 
196
    return normalize(b, a)
 
197
 
 
198
def lp2hp(b,a,wo=1.0):
 
199
    """Return a high-pass filter with cuttoff frequency wo
 
200
    from a low-pass filter prototype with unity cutoff frequency.
 
201
    """
 
202
    a,b = map(atleast_1d,(a,b))
 
203
    if type(wo) is type(a):
 
204
        wo = wo[0]    
 
205
    d = len(a)
 
206
    n = len(b)
 
207
    if wo != 1:
 
208
        pwo = pow(wo,Num.arange(max((d,n))))
 
209
    else:
 
210
        pwo = Num.ones(max((d,n)),b.typecode())
 
211
    if d >= n:
 
212
        outa = a[::-1] * pwo
 
213
        outb = resize(b,(d,))
 
214
        outb[n:] = 0.0
 
215
        outb[:n] = b[::-1] * pwo[:n]        
 
216
    else:
 
217
        outb = b[::-1] * pwo
 
218
        outa = resize(a,(n,))
 
219
        outa[d:] = 0.0
 
220
        outa[:d] = a[::-1] * pwo[:d]
 
221
 
 
222
    return normalize(outb, outa)
 
223
 
 
224
def lp2bp(b,a,wo=1.0, bw=1.0):
 
225
    """Return a band-pass filter with center frequency wo and bandwidth bw
 
226
    from a low-pass filter prototype with unity cutoff frequency.
 
227
    """
 
228
    a,b = map(atleast_1d,(a,b))
 
229
    D = len(a) - 1
 
230
    N = len(b) - 1
 
231
    artype = b.typecode()
 
232
    if artype not in ['F','D','f','d']:
 
233
        artype = Num.Float
 
234
    ma = max([N,D])
 
235
    Np = N + ma
 
236
    Dp = D + ma
 
237
    bprime = Num.zeros(Np+1,artype)
 
238
    aprime = Num.zeros(Dp+1,artype)
 
239
    wosq = wo*wo
 
240
    for j in range(Np+1):
 
241
        val = 0.0
 
242
        for i in range(0,N+1):
 
243
            for k in range(0,i+1):
 
244
                if ma-i+2*k == j:
 
245
                    val += comb(i,k)*b[N-i]*(wosq)**(i-k) / bw**i
 
246
        bprime[Np-j] = val
 
247
    for j in range(Dp+1):
 
248
        val = 0.0
 
249
        for i in range(0,D+1):
 
250
            for k in range(0,i+1):
 
251
                if ma-i+2*k == j:
 
252
                    val += comb(i,k)*a[D-i]*(wosq)**(i-k) / bw**i
 
253
        aprime[Dp-j] = val
 
254
        
 
255
    return normalize(bprime, aprime)
 
256
 
 
257
def lp2bs(b,a,wo=1,bw=1):
 
258
    """Return a band-stop filter with center frequency wo and bandwidth bw
 
259
    from a low-pass filter prototype with unity cutoff frequency.
 
260
    """
 
261
    a,b = map(atleast_1d,(a,b))
 
262
    D = len(a) - 1
 
263
    N = len(b) - 1
 
264
    artype = b.typecode()
 
265
    if artype not in ['F','D','f','d']:
 
266
        artype = Num.Float
 
267
    M = max([N,D])
 
268
    Np = M + M
 
269
    Dp = M + M
 
270
    bprime = Num.zeros(Np+1,artype)
 
271
    aprime = Num.zeros(Dp+1,artype)
 
272
    wosq = wo*wo
 
273
    for j in range(Np+1):
 
274
        val = 0.0
 
275
        for i in range(0,N+1):
 
276
            for k in range(0,M-i+1):
 
277
                if i+2*k == j:
 
278
                    val += comb(M-i,k)*b[N-i]*(wosq)**(M-i-k) * bw**i
 
279
        bprime[Np-j] = val
 
280
    for j in range(Dp+1):
 
281
        val = 0.0
 
282
        for i in range(0,D+1):
 
283
            for k in range(0,M-i+1):
 
284
                if i+2*k == j:
 
285
                    val += comb(M-i,k)*a[D-i]*(wosq)**(M-i-k) * bw**i
 
286
        aprime[Dp-j] = val
 
287
        
 
288
    return normalize(bprime, aprime)
 
289
 
 
290
def bilinear(b,a,fs=1.0):
 
291
    """Return a digital filter from an analog filter using the bilinear transform.
 
292
 
 
293
    The bilinear transform substitutes (z-1) / (z+1) for s
 
294
    """
 
295
    fs =float(fs)
 
296
    a,b = map(atleast_1d,(a,b))
 
297
    D = len(a) - 1
 
298
    N = len(b) - 1
 
299
    artype = Num.Float
 
300
    M = max([N,D])
 
301
    Np = M
 
302
    Dp = M
 
303
    bprime = Num.zeros(Np+1,artype)
 
304
    aprime = Num.zeros(Dp+1,artype)
 
305
    for j in range(Np+1):
 
306
        val = 0.0
 
307
        for i in range(N+1):
 
308
             for k in range(i+1):
 
309
                for l in range(M-i+1):
 
310
                    if k+l == j:
 
311
                        val += comb(i,k)*comb(M-i,l)*b[N-i]*pow(2*fs,i)*(-1)**k
 
312
        bprime[j] = real(val)
 
313
    for j in range(Dp+1):
 
314
        val = 0.0
 
315
        for i in range(D+1):
 
316
            for k in range(i+1):
 
317
                for l in range(M-i+1):
 
318
                    if k+l == j:
 
319
                        val += comb(i,k)*comb(M-i,l)*a[D-i]*pow(2*fs,i)*(-1)**k
 
320
        aprime[j] = real(val)
 
321
        
 
322
    return normalize(bprime, aprime)
 
323
 
 
324
def iirdesign(wp, ws, gpass, gstop, analog=0, ftype='ellip', output='ba'):
 
325
    """Complete IIR digital and analog filter design.
 
326
 
 
327
    Description:
 
328
 
 
329
      Given passband and stopband frequencies and gains construct an analog or
 
330
      digital IIR filter of minimum order for a given basic type.  Return the
 
331
      output in numerator, denominator ('ba') or pole-zero ('zpk') form.
 
332
 
 
333
    Inputs:
 
334
 
 
335
      wp, ws -- Passband and stopband edge frequencies, normalized from 0
 
336
                to 1 (1 corresponds to pi radians / sample).  For example:
 
337
                   Lowpass:   wp = 0.2,          ws = 0.3
 
338
                   Highpass:  wp = 0.3,          ws = 0.2
 
339
                   Bandpass:  wp = [0.2, 0.5],   ws = [0.1, 0.6]
 
340
                   Bandstop:  wp = [0.1, 0.6],   ws = [0.2, 0.5]
 
341
      gpass -- The maximum loss in the passband (dB).
 
342
      gstop -- The minimum attenuation in the stopband (dB).
 
343
      analog -- Non-zero to design an analog filter (in this case wp and
 
344
                ws are in radians / second).
 
345
      ftype -- The type of iir filter to design:
 
346
                 elliptic    : 'ellip'
 
347
                 Butterworth : 'butter',
 
348
                 Chebyshev I : 'cheby1',
 
349
                 Chebyshev II: 'cheby2',
 
350
                 Bessel :      'bessel'
 
351
      output -- Type of output:  numerator/denominator ('ba') or pole-zero ('zpk')
 
352
 
 
353
    Outputs: (b,a) or (z,p,k)
 
354
 
 
355
      b,a -- Numerator and denominator of the iir filter.
 
356
      z,p,k -- Zeros, poles, and gain of the iir filter.      
 
357
    """
 
358
 
 
359
    try:
 
360
        ordfunc = filter_dict[ftype][1]
 
361
    except KeyError:
 
362
        raise ValueError, "Invalid IIR filter type."
 
363
    except IndexError:
 
364
        raise ValueError, "%s does not have order selection use iirfilter function." % ftype
 
365
 
 
366
    wp = atleast_1d(wp)
 
367
    ws = atleast_1d(ws)
 
368
    band_type = 2*(len(wp)-1)
 
369
    band_type +=1 
 
370
    if wp[0] >= ws[0]:
 
371
        band_type += 1
 
372
 
 
373
    btype = {1:'lowpass', 2:'highpass', 3:'bandstop', 4:'bandpass'}[band_type]
 
374
       
 
375
    N, Wn = ordfunc(wp, ws, gpass, gstop, analog=analog)
 
376
    return iirfilter(N, Wn, rp=gpass, rs=gstop, analog=analog, btype=btype, ftype=ftype, output=output)
 
377
 
 
378
 
 
379
def iirfilter(N, Wn, rp=None, rs=None, btype='band', analog=0, ftype='butter', output='ba'):
 
380
    """IIR digital and analog filter design given order and critical points.
 
381
 
 
382
    Description:
 
383
 
 
384
      Design an Nth order lowpass digital or analog filter and return the filter
 
385
      coefficients in (B,A) (numerator, denominator) or (Z,P,K) form.          
 
386
 
 
387
    Inputs:
 
388
 
 
389
      N -- the order of the filter.
 
390
      Wn -- a scalar or length-2 sequence giving the critical frequencies.
 
391
      rp, rs -- For chebyshev and elliptic filters provides the maximum ripple
 
392
                in the passband and the minimum attenuation in the stop band.
 
393
      btype -- the type of filter (lowpass, highpass, bandpass, or bandstop).
 
394
      analog -- non-zero to return an analog filter, otherwise
 
395
                a digital filter is returned.
 
396
      ftype -- the type of IIR filter (Butterworth, Cauer (Elliptic),
 
397
               Bessel, Chebyshev1, Chebyshev2)
 
398
      output -- 'ba' for (b,a) output, 'zpk' for (z,p,k) output.
 
399
 
 
400
    SEE ALSO butterord, cheb1ord, cheb2ord, ellipord
 
401
    """
 
402
 
 
403
    ftype, btype, output = map(string.lower, (ftype, btype, output))
 
404
    Wn = asarray(Wn)
 
405
    try:
 
406
        btype = band_dict[btype]
 
407
    except KeyError:
 
408
        raise ValueError, "%s is an invalid bandtype for filter." % btype
 
409
 
 
410
    try:
 
411
        typefunc = filter_dict[ftype][0]
 
412
    except KeyError:
 
413
        raise ValueError, "%s is not a valid basic iir filter." % ftype
 
414
 
 
415
    if output not in ['ba', 'zpk']:
 
416
        raise ValueError, "%s is not a valid output form." % output
 
417
    
 
418
    #pre-warp frequencies for digital filter design
 
419
    if not analog:
 
420
        fs = 2.0
 
421
        warped = 2*fs*tan(pi*Wn/fs)
 
422
    else:
 
423
        warped = Wn
 
424
    
 
425
    # convert to low-pass prototype
 
426
    if btype in ['lowpass', 'highpass']:
 
427
        wo = warped
 
428
    else:
 
429
        bw = warped[1] - warped[0]
 
430
        wo = sqrt(warped[0]*warped[1])    
 
431
    
 
432
    # Get analog lowpass prototype
 
433
    if typefunc in [buttap, besselap]:
 
434
        z, p, k = typefunc(N)
 
435
    elif typefunc == cheb1ap:
 
436
        if rp is None:
 
437
            raise ValueError, "passband ripple (rp) must be provided to design a Chebyshev I filter."
 
438
        z, p, k = typefunc(N, rp)
 
439
    elif typefunc == cheb2ap:
 
440
        if rs is None:
 
441
            raise ValueError, "stopband atteunatuion (rs) must be provided to design an Chebyshev II filter."
 
442
        z, p, k = typefunc(N, rs)
 
443
    else:  # Elliptic filters
 
444
        if rs is None or rp is None:
 
445
            raise ValueErrro, "Both rp and rs must be provided to design an elliptic filter."
 
446
        z, p, k = typefunc(N, rp, rs)
 
447
    
 
448
    b, a = zpk2tf(z,p,k)
 
449
    
 
450
    # transform to lowpass, bandpass, highpass, or bandstop
 
451
    if btype == 'lowpass':
 
452
        b, a = lp2lp(b,a,wo=wo)
 
453
    elif btype == 'highpass':
 
454
        b, a = lp2hp(b,a,wo=wo)
 
455
    elif btype == 'bandpass':
 
456
        b, a = lp2bp(b,a,wo=wo,bw=bw)
 
457
    else: # 'bandstop'
 
458
        b, a = lp2bs(b,a,wo=wo,bw=bw)
 
459
               
 
460
 
 
461
    # Find discrete equivalent if necessary
 
462
    if not analog:
 
463
        b, a = bilinear(b, a, fs=fs)
 
464
    
 
465
    # Transform to proper out type (pole-zero, state-space, numer-denom)    
 
466
    if output == 'zpk':
 
467
        return tf2zpk(b,a)
 
468
    else:
 
469
        return b,a
 
470
      
 
471
 
 
472
def butter(N, Wn, btype='low', analog=0, output='ba'):
 
473
    """Butterworth digital and analog filter design.
 
474
 
 
475
    Description:
 
476
 
 
477
      Design an Nth order lowpass digital or analog Butterworth filter
 
478
      and return the filter coefficients in (B,A) or (Z,P,K) form.
 
479
 
 
480
    See also buttord.
 
481
    """
 
482
    return iirfilter(N, Wn, btype=btype, analog=analog, output=output, ftype='butter')
 
483
 
 
484
def cheby1(N, rp, Wn, btype='low', analog=0, output='ba'):
 
485
    """Chebyshev type I digital and analog filter design.
 
486
 
 
487
    Description:
 
488
 
 
489
      Design an Nth order lowpass digital or analog Chebyshev type I filter
 
490
      and return the filter coefficients in (B,A) or (Z,P,K) form.
 
491
 
 
492
    See also cheb1ord.      
 
493
    """
 
494
    return iirfilter(N, Wn, rp=rp, btype=btype, analog=analog, output=output, ftype='cheby1')
 
495
 
 
496
def cheby2(N, rs, Wn, btype='low', analog=0, output='ba'):
 
497
    """Chebyshev type I digital and analog filter design.
 
498
 
 
499
    Description:
 
500
 
 
501
      Design an Nth order lowpass digital or analog Chebyshev type I filter
 
502
      and return the filter coefficients in (B,A) or (Z,P,K) form.
 
503
 
 
504
    See also cheb2ord.
 
505
    """
 
506
    return iirfilter(N, Wn, rs=rs, btype=btype, analog=analog, output=output, ftype='cheby2')
 
507
 
 
508
def ellip(N, rp, rs, Wn, btype='low', analog=0, output='ba'):
 
509
    """Elliptic (Cauer) digital and analog filter design.
 
510
 
 
511
    Description:
 
512
 
 
513
      Design an Nth order lowpass digital or analog elliptic filter
 
514
      and return the filter coefficients in (B,A) or (Z,P,K) form.
 
515
 
 
516
    See also ellipord.
 
517
    """
 
518
    return iirfilter(N, Wn, rs=rs, rp=rp, btype=btype, analog=analog, output=output, ftype='elliptic')
 
519
 
 
520
def bessel(N, Wn, btype='low', analog=0, output='ba'):
 
521
    """Bessel digital and analog filter design.
 
522
 
 
523
    Description:
 
524
 
 
525
      Design an Nth order lowpass digital or analog Bessel filter
 
526
      and return the filter coefficients in (B,A) or (Z,P,K) form.
 
527
 
 
528
    """
 
529
    return iirfilter(N, Wn, btype=btype, analog=analog, output=output, ftype='bessel')
 
530
 
 
531
    
 
532
def maxflat():
 
533
    pass
 
534
 
 
535
def yulewalk():
 
536
    pass
 
537
 
 
538
 
 
539
def band_stop_obj(wp, ind, passb, stopb, gpass, gstop, type):
 
540
    """Band Stop Objective Function for order minimization
 
541
 
 
542
    Description:
 
543
 
 
544
      Returns the non-integer order for an analog band stop filter.
 
545
 
 
546
    Inputs:
 
547
 
 
548
      wp -- passb edge
 
549
      ind -- index specifying which passb edge to vary (0 or 1).
 
550
      passb -- two element vector of fixed passband edges.
 
551
      stopb -- two element vector of fixed stopband edges.
 
552
      gstop -- amount in dB of attenuation in stopband.
 
553
      gpass -- amount in dB of ripple in the passband.
 
554
      type -- 'butter', 'cheby', or 'ellip':
 
555
 
 
556
    Outputs: (n,)
 
557
 
 
558
      n -- filter order (possibly non-integer)
 
559
    """
 
560
 
 
561
    passbC = passb.copy()
 
562
    passbC[ind] = wp
 
563
    nat = stopb*(passbC[0]-passbC[1]) / (stopb**2 - passbC[0]*passbC[1])
 
564
    nat = min(abs(nat))
 
565
 
 
566
    if type == 'butter':
 
567
        GSTOP = 10**(0.1*abs(gstop))
 
568
        GPASS = 10**(0.1*abs(gpass))
 
569
        n = (log10((GSTOP-1.0)/(GPASS-1.0)) / (2*log10(nat)))
 
570
    elif type == 'cheby':
 
571
        GSTOP = 10**(0.1*abs(gstop))
 
572
        GPASS = 10**(0.1*abs(gpass))
 
573
        n = arccosh(sqrt((GSTOP-1.0)/(GPASS-1.0))) / arccosh(nat)
 
574
    elif type == 'ellip':
 
575
        GSTOP = 10**(0.1*gstop)
 
576
        GPASS = 10**(0.1*gpass)
 
577
        arg1 = sqrt( (GPASS-1.0) / (GSTOP-1.0) )
 
578
        arg0 = 1.0 / nat
 
579
        d0 = special.ellipk([arg0**2, 1-arg0**2])
 
580
        d1 = special.ellipk([arg1**2, 1-arg1**2])
 
581
        n = (d0[0]*d1[1] / (d0[1]*d1[0]))
 
582
    else:
 
583
        raise ValueError, "Incorrect type: ", type
 
584
    return n
 
585
 
 
586
def buttord(wp, ws, gpass, gstop, analog=0):
 
587
    """Butterworth filter order selection.
 
588
 
 
589
    Description:
 
590
 
 
591
      Return the order of the lowest order digital Butterworth filter that
 
592
      loses no more than gpass dB in the passband and has at least gstop dB
 
593
      attenuation in the stopband.
 
594
 
 
595
    Inputs:
 
596
 
 
597
      wp, ws -- Passband and stopband edge frequencies, normalized from 0
 
598
                to 1 (1 corresponds to pi radians / sample).  For example:
 
599
                   Lowpass:   wp = 0.2,          ws = 0.3
 
600
                   Highpass:  wp = 0.3,          ws = 0.2
 
601
                   Bandpass:  wp = [0.2, 0.5],   ws = [0.1, 0.6]
 
602
                   Bandstop:  wp = [0.1, 0.6],   ws = [0.2, 0.5]
 
603
      gpass -- The maximum loss in the passband (dB).
 
604
      gstop -- The minimum attenuation in the stopband (dB).
 
605
      analog -- Non-zero to design an analog filter (in this case wp and
 
606
                ws are in radians / second).
 
607
 
 
608
    Outputs: (ord, Wn)
 
609
 
 
610
      ord -- The lowest order for a Butterworth filter which meets specs.
 
611
      Wn -- The Butterworth natural frequency (i.e. the "3dB frequency"). 
 
612
            Should be used with scipy.signal.butter to give filter results.
 
613
 
 
614
    """
 
615
 
 
616
    wp = atleast_1d(wp)
 
617
    ws = atleast_1d(ws)
 
618
    filter_type = 2*(len(wp)-1)
 
619
    filter_type +=1 
 
620
    if wp[0] >= ws[0]:
 
621
        filter_type += 1
 
622
 
 
623
    # Pre-warp frequencies
 
624
    if not analog:
 
625
        passb = tan(wp*pi/2.0)
 
626
        stopb = tan(ws*pi/2.0)
 
627
    else:
 
628
        passb = wp
 
629
        stopb = ws
 
630
 
 
631
    if filter_type == 1:            # low
 
632
        nat = stopb / passb
 
633
    elif filter_type == 2:          # high
 
634
        nat = passb / stopb
 
635
    elif filter_type == 3:          # stop
 
636
        wp0 = optimize.fminbound(band_stop_obj, passb[0], stopb[0]-1e-12,
 
637
                                 args=(0,passb,stopb,gpass,gstop,'butter'),
 
638
                                 disp=0)
 
639
        passb[0] = wp0
 
640
        wp1 = optimize.fminbound(band_stop_obj, stopb[1]+1e-12, passb[1],
 
641
                                 args=(1,passb,stopb,gpass,gstop,'butter'),
 
642
                                 disp=0)
 
643
        passb[1] = wp1
 
644
        nat = (stopb * (passb[0]-passb[1])) / (stopb**2 - passb[0]*passb[1])
 
645
    elif filter_type == 4:          # pass
 
646
        nat = (stopb**2 - passb[0]*passb[1]) / (stopb* (passb[0]-passb[1]))
 
647
 
 
648
    nat = min(abs(nat))
 
649
 
 
650
    GSTOP = 10**(0.1*abs(gstop))
 
651
    GPASS = 10**(0.1*abs(gpass))
 
652
    ord = int(ceil( log10((GSTOP-1.0)/(GPASS-1.0)) / (2*log10(nat))))
 
653
 
 
654
    # Find the butterworth natural frequency W0 (or the "3dB" frequency")
 
655
    # to give exactly gstop at nat. W0 will be between 1 and nat
 
656
    try:
 
657
        W0 = nat / ( ( 10**(0.1*abs(gstop))-1)**(1.0/(2.0*ord)))
 
658
    except ZeroDivisionError:
 
659
        W0 = nat
 
660
        print "Warning, order is zero...check input parametegstop."
 
661
 
 
662
    # now convert this frequency back from lowpass prototype
 
663
    # to the original analog filter
 
664
 
 
665
    if filter_type == 1:  # low
 
666
        WN = W0*passb
 
667
    elif filter_type == 2: # high
 
668
        WN = passb / W0
 
669
    elif filter_type == 3:  # stop
 
670
        WN = Num.zeros(2,Float)
 
671
        WN[0] = ((passb[1] - passb[0]) + sqrt((passb[1] - passb[0])**2 + \
 
672
                                        4*W0**2 * passb[0] * passb[1])) / (2*W0)
 
673
        WN[1] = ((passb[1] - passb[0]) - sqrt((passb[1] - passb[0])**2 + \
 
674
                                        4*W0**2 * passb[0] * passb[1])) / (2*W0)
 
675
        WN = Num.sort(abs(WN))
 
676
    elif filter_type == 4: # pass
 
677
        W0 = Num.array([-W0, W0],Float)
 
678
        WN = -W0 * (passb[1]-passb[0]) / 2.0 + sqrt(W0**2 / 4.0 * \
 
679
                                              (passb[1]-passb[0])**2 + \
 
680
                                              passb[0]*passb[1])
 
681
        WN = Num.sort(abs(WN))
 
682
    else:
 
683
        raise ValueError, "Bad type."
 
684
 
 
685
    if not analog:
 
686
        wn = (2.0/pi)*arctan(WN)
 
687
    else:
 
688
        wn = WN
 
689
 
 
690
    if len(wn) == 1:
 
691
        wn = wn[0]
 
692
    return ord, wn
 
693
 
 
694
def cheb1ord(wp, ws, gpass, gstop, analog=0):
 
695
    """Chebyshev type I filter order selection.
 
696
 
 
697
    Description:
 
698
 
 
699
      Return the order of the lowest order digital Chebyshev Type I filter
 
700
      that loses no more than gpass dB in the passband and has at least gstop dB
 
701
      attenuation in the stopband.
 
702
 
 
703
    Inputs:
 
704
 
 
705
      wp, ws -- Passband and stopband edge frequencies, normalized from 0
 
706
                to 1 (1 corresponds to pi radians / sample).  For example:
 
707
                   Lowpass:   wp = 0.2,          ws = 0.3
 
708
                   Highpass:  wp = 0.3,          ws = 0.2
 
709
                   Bandpass:  wp = [0.2, 0.5],   ws = [0.1, 0.6]
 
710
                   Bandstop:  wp = [0.1, 0.6],   ws = [0.2, 0.5]
 
711
      gpass -- The maximum loss in the passband (dB).
 
712
      gstop -- The minimum attenuation in the stopband (dB).
 
713
      analog -- Non-zero to design an analog filter (in this case wp and
 
714
                ws are in radians / second).
 
715
 
 
716
    Outputs: (ord, Wn)
 
717
 
 
718
      ord -- The lowest order for a Chebyshev type I filter that meets specs.
 
719
      Wn -- The Chebyshev natural frequency (the "3dB frequency") for
 
720
            use with scipy.signal.cheby1 to give filter results.
 
721
 
 
722
    """
 
723
    wp = atleast_1d(wp)
 
724
    ws = atleast_1d(ws)
 
725
    filter_type = 2*(len(wp)-1)
 
726
    if wp[0] < ws[0]:
 
727
        filter_type += 1
 
728
    else:
 
729
        filter_type += 2
 
730
 
 
731
    # Pre-wagpass frequencies
 
732
    if not analog:
 
733
        passb = tan(pi*wp/2)
 
734
        stopb = tan(pi*ws/2)
 
735
    else:
 
736
        passb = wp
 
737
        stopb = ws
 
738
 
 
739
    if filter_type == 1:           # low
 
740
        nat = stopb / passb
 
741
    elif filter_type == 2:          # high
 
742
        nat = passb / stopb
 
743
    elif filter_type == 3:     # stop
 
744
        wp0 = optimize.fminbound(band_stop_obj, passb[0], stopb[0]-1e-12,
 
745
                                 args=(0,passb,stopb,gpass,gstop,'cheby'), disp=0)
 
746
        passb[0] = wp0
 
747
        wp1 = optimize.fminbound(band_stop_obj, stopb[1]+1e-12, passb[1],
 
748
                                 args=(1,passb,stopb,gpass,gstop,'cheby'), disp=0)
 
749
        passb[1] = wp1
 
750
        nat = (stopb * (passb[0]-passb[1])) / (stopb**2 - passb[0]*passb[1])
 
751
    elif filter_type == 4:  # pass
 
752
        nat = (stopb**2 - passb[0]*passb[1]) / (stopb* (passb[0]-passb[1]))
 
753
 
 
754
    nat = min(abs(nat))
 
755
 
 
756
    GSTOP = 10**(0.1*abs(gstop))
 
757
    GPASS = 10**(0.1*abs(gpass))
 
758
    ord = int(ceil(arccosh(sqrt((GSTOP-1.0) / (GPASS-1.0))) / arccosh(nat)))
 
759
 
 
760
    # Natural frequencies are just the passband edges 
 
761
    if not analog:
 
762
        wn = (2.0/pi)*arctan(passb)
 
763
    else:
 
764
        wn = passb
 
765
 
 
766
    if len(wn) == 1:
 
767
        wn = wn[0]
 
768
    return ord, wn
 
769
    
 
770
 
 
771
def cheb2ord(wp, ws, gpass, gstop, analog=0):
 
772
    """Chebyshev type II filter order selection.
 
773
 
 
774
    Description:
 
775
 
 
776
      Return the order of the lowest order digital Chebyshev Type II filter
 
777
      that loses no more than gpass dB in the passband and has at least gstop dB
 
778
      attenuation in the stopband.
 
779
 
 
780
    Inputs:
 
781
 
 
782
      wp, ws -- Passband and stopband edge frequencies, normalized from 0
 
783
                to 1 (1 corresponds to pi radians / sample).  For example:
 
784
                   Lowpass:   wp = 0.2,          ws = 0.3
 
785
                   Highpass:  wp = 0.3,          ws = 0.2
 
786
                   Bandpass:  wp = [0.2, 0.5],   ws = [0.1, 0.6]
 
787
                   Bandstop:  wp = [0.1, 0.6],   ws = [0.2, 0.5]
 
788
      gpass -- The maximum loss in the passband (dB).
 
789
      gstop -- The minimum attenuation in the stopband (dB).
 
790
      analog -- Non-zero to design an analog filter (in this case wp and
 
791
                ws are in radians / second).
 
792
 
 
793
    Outputs: (ord, Wn)
 
794
 
 
795
      ord -- The lowest order for a Chebyshev type II filter that meets specs.
 
796
      Wn -- The Chebyshev natural frequency for
 
797
            use with scipy.signal.cheby2 to give the filter.
 
798
 
 
799
    """
 
800
    wp = atleast_1d(wp)
 
801
    ws = atleast_1d(ws)
 
802
    filter_type = 2*(len(wp)-1)
 
803
    if wp[0] < ws[0]:
 
804
        filter_type += 1
 
805
    else:
 
806
        filter_type += 2
 
807
 
 
808
    # Pre-wagpass frequencies
 
809
    if not analog:
 
810
        passb = tan(pi*wp/2)
 
811
        stopb = tan(pi*ws/2)
 
812
    else:
 
813
        passb = wp
 
814
        stopb = ws
 
815
 
 
816
    if filter_type == 1:           # low
 
817
        nat = stopb / passb
 
818
    elif filter_type == 2:          # high
 
819
        nat = passb / stopb
 
820
    elif filter_type == 3:     # stop
 
821
        wp0 = optimize.fminbound(band_stop_obj, passb[0], stopb[0]-1e-12,
 
822
                                 args=(0,passb,stopb,gpass,gstop,'cheby'),
 
823
                                 disp=0)
 
824
        passb[0] = wp0
 
825
        wp1 = optimize.fminbound(band_stop_obj, stopb[1]+1e-12, passb[1],
 
826
                                 args=(1,passb,stopb,gpass,gstop,'cheby'),
 
827
                                 disp=0)
 
828
        passb[1] = wp1
 
829
        nat = (stopb * (passb[0]-passb[1])) / (stopb**2 - passb[0]*passb[1])
 
830
    elif filter_type == 4:  # pass
 
831
        nat = (stopb**2 - passb[0]*passb[1]) / (stopb* (passb[0]-passb[1]))
 
832
 
 
833
    nat = min(abs(nat))
 
834
 
 
835
    GSTOP = 10**(0.1*abs(gstop))
 
836
    GPASS = 10**(0.1*abs(gpass))
 
837
    ord = int(ceil(arccosh(sqrt((GSTOP-1.0) / (GPASS-1.0))) / arccosh(nat)))
 
838
 
 
839
    # Find frequency where analog response is -gpass dB.
 
840
    # Then convert back from low-pass prototype to the original filter.
 
841
 
 
842
    new_freq = cosh(1.0/ord * arccosh(sqrt((GSTOP-1.0)/(GPASS-1.0))))
 
843
    new_freq = 1.0 / new_freq
 
844
    
 
845
    if filter_type == 1:
 
846
        nat = passb / new_freq
 
847
    elif filter_type == 2:
 
848
        nat = passb * new_freq
 
849
    elif filter_type == 3:
 
850
        nat = Num.zeros(2,Num.Float)
 
851
        nat[0] = new_freq / 2.0 * (passb[0]-passb[1]) + \
 
852
                 sqrt(new_freq**2 * (passb[1]-passb[0])**2 / 4.0 + \
 
853
                      passb[1] * passb[0])
 
854
        nat[1] = passb[1] * passb[0] / nat[0]
 
855
    elif filter_type == 4:
 
856
        nat = Num.zeros(2,Num.Float)
 
857
        nat[0] = 1.0/(2.0*new_freq) * (passb[0] - passb[1]) + \
 
858
                 sqrt((passb[1]-passb[0])**2 / (4.0*new_freq**2) + \
 
859
                      passb[1] * passb[0])
 
860
        nat[1] = passb[0] * passb[1] / nat[0]        
 
861
 
 
862
    if not analog:
 
863
        wn = (2.0/pi)*arctan(nat)
 
864
    else:
 
865
        wn = nat
 
866
 
 
867
    if len(wn) == 1:
 
868
        wn = wn[0]
 
869
    return ord, wn
 
870
 
 
871
 
 
872
def ellipord(wp, ws, gpass, gstop, analog=0):
 
873
    """Elliptic (Cauer) filter order selection.
 
874
 
 
875
    Description:
 
876
 
 
877
      Return the order of the lowest order digital elliptic filter
 
878
      that loses no more than gpass dB in the passband and has at least gstop dB
 
879
      attenuation in the stopband.
 
880
 
 
881
    Inputs:
 
882
 
 
883
      wp, ws -- Passband and stopband edge frequencies, normalized from 0
 
884
                to 1 (1 corresponds to pi radians / sample).  For example:
 
885
                   Lowpass:   wp = 0.2,          ws = 0.3
 
886
                   Highpass:  wp = 0.3,          ws = 0.2
 
887
                   Bandpass:  wp = [0.2, 0.5],   ws = [0.1, 0.6]
 
888
                   Bandstop:  wp = [0.1, 0.6],   ws = [0.2, 0.5]
 
889
      gpass -- The maximum loss in the passband (dB).
 
890
      gstop -- The minimum attenuation in the stopband (dB).
 
891
      analog -- Non-zero to design an analog filter (in this case wp and
 
892
                ws are in radians / second).
 
893
 
 
894
    Outputs: (ord, Wn)
 
895
 
 
896
      ord -- The lowest order for a Chebyshev type II filter that meets specs.
 
897
      Wn -- The Chebyshev natural frequency for
 
898
            use with scipy.signal.cheby2 to give the filter.
 
899
 
 
900
    """
 
901
    wp = atleast_1d(wp)
 
902
    ws = atleast_1d(ws)
 
903
    filter_type = 2*(len(wp)-1)
 
904
    filter_type += 1
 
905
    if wp[0] >= ws[0]:
 
906
        filter_type += 1
 
907
 
 
908
    # Pre-wagpass frequencies
 
909
    if analog:
 
910
        passb = wp
 
911
        stopb = ws
 
912
    else:
 
913
        passb = tan(wp*pi/2.0)
 
914
        stopb = tan(ws*pi/2.0)
 
915
 
 
916
    if filter_type == 1:           # low
 
917
        nat = stopb / passb
 
918
    elif filter_type == 2:          # high
 
919
        nat = passb / stopb
 
920
    elif filter_type == 3:     # stop
 
921
        wp0 = optimize.fminbound(band_stop_obj, passb[0], stopb[0]-1e-12,
 
922
                                 args=(0,passb,stopb,gpass,gstop,'ellip'),
 
923
                                 disp=0)
 
924
        passb[0] = wp0
 
925
        wp1 = optimize.fminbound(band_stop_obj, stopb[1]+1e-12, passb[1],
 
926
                                 args=(1,passb,stopb,gpass,gstop,'ellip'),
 
927
                                 disp=0)
 
928
        passb[1] = wp1
 
929
        nat = (stopb * (passb[0]-passb[1])) / (stopb**2 - passb[0]*passb[1])
 
930
    elif filter_type == 4:  # pass
 
931
        nat = (stopb**2 - passb[0]*passb[1]) / (stopb* (passb[0]-passb[1]))
 
932
 
 
933
    nat = min(abs(nat))
 
934
 
 
935
    GSTOP = 10**(0.1*gstop)
 
936
    GPASS = 10**(0.1*gpass)
 
937
    arg1 = sqrt( (GPASS-1.0) / (GSTOP-1.0) )
 
938
    arg0 = 1.0 / nat
 
939
    d0 = special.ellipk([arg0**2, 1-arg0**2])
 
940
    d1 = special.ellipk([arg1**2, 1-arg1**2])
 
941
    ord = int(ceil(d0[0]*d1[1] / (d0[1]*d1[0])))
 
942
 
 
943
    if not analog:
 
944
        wn = arctan(passb)*2.0/pi
 
945
    else:
 
946
        wn = passb
 
947
 
 
948
    if len(wn) == 1:
 
949
        wn = wn[0]
 
950
    return ord, wn
 
951
    
 
952
def buttap(N):
 
953
    """Return (z,p,k) zero, pole, gain for analog prototype of an Nth
 
954
    order Butterworth filter."""
 
955
    z = []
 
956
    n = Num.arange(1,N+1)
 
957
    p = Num.exp(1j*(2*n-1)/(2.0*N)*pi)*1j
 
958
    k = 1
 
959
    return z, p, k
 
960
 
 
961
def cheb1ap(N,rp):
 
962
    """Return (z,p,k) zero, pole, gain for Nth order Chebyshev type I
 
963
    lowpass analog filter prototype with rp decibels of ripple
 
964
    in the passband.
 
965
    """
 
966
    z = []
 
967
    eps = Num.sqrt(10**(0.1*rp)-1.0)
 
968
    n = Num.arange(1,N+1)
 
969
    mu = 1.0/N * Num.log((1.0+Num.sqrt(1+eps*eps)) / eps)
 
970
    theta = pi/2.0 * (2*n-1.0)/N
 
971
    p = -Num.sinh(mu)*Num.sin(theta) + 1j*Num.cosh(mu)*Num.cos(theta)
 
972
    k = MLab.prod(-p).real
 
973
    if N % 2 == 0:
 
974
        k = k / sqrt((1+eps*eps))
 
975
    return z, p, k
 
976
    pass
 
977
 
 
978
def cheb2ap(N,rs):
 
979
    """Return (z,p,k) zero, pole, gain for Nth order Chebyshev type II
 
980
    lowpass analog filter prototype with rs decibels of ripple
 
981
    in the stopband.
 
982
    """
 
983
    de = 1.0/sqrt(10**(0.1*rs)-1)
 
984
    mu = arcsinh(1.0/de)/N
 
985
 
 
986
    if N % 2:
 
987
        m = N - 1
 
988
        n = Num.concatenate((Num.arange(1,N-1,2),Num.arange(N+2,2*N,2)))
 
989
    else:
 
990
        m = N
 
991
        n = Num.arange(1,2*N,2)
 
992
        
 
993
    z = conjugate(1j / cos(n*pi/(2.0*N)))
 
994
    p = exp(1j*(pi*Num.arange(1,2*N,2)/(2.0*N) + pi/2.0))
 
995
    p = sinh(mu) * p.real + 1j*cosh(mu)*p.imag
 
996
    p = 1.0 / p
 
997
    k = (MLab.prod(-p)/MLab.prod(-z)).real
 
998
    return z, p, k
 
999
    
 
1000
 
 
1001
EPSILON = 2e-16
 
1002
 
 
1003
def vratio(u, ineps, mp):
 
1004
    [s,c,d,phi] = special.ellipj(u,mp)
 
1005
    ret = abs(ineps - s/c)
 
1006
    return ret
 
1007
 
 
1008
def kratio(m, k_ratio):
 
1009
    m = float(m)
 
1010
    if m < 0:
 
1011
        m = 0.0
 
1012
    if m > 1:
 
1013
        m = 1.0
 
1014
    if abs(m) > EPSILON and (abs(m) + EPSILON) < 1:
 
1015
        k = special.ellipk([m,1-m])
 
1016
        r = k[0] / k[1] - k_ratio
 
1017
    elif abs(m) > EPSILON:
 
1018
        r = -k_ratio
 
1019
    else:
 
1020
        r = 1e20
 
1021
    return abs(r)
 
1022
 
 
1023
def ellipap(N,rp,rs):
 
1024
    """Return (z,p,k) zeros, poles, and gain of an Nth order normalized
 
1025
    prototype elliptic analog lowpass filter with rp decibels of ripple
 
1026
    in the passband and a stopband rs decibels down.
 
1027
 
 
1028
    See Chapter 12 and Chapter 5 of "Filter Design for Signal Processing",
 
1029
    by Lutova, Tosic, and Evans.  This is 
 
1030
    """
 
1031
    if N == 1:
 
1032
        p = -sqrt(1.0/(10**(0.1*rp)-1.0))
 
1033
        k = -p
 
1034
        z = []
 
1035
        return z, p, k
 
1036
 
 
1037
    eps = Num.sqrt(10**(0.1*rp)-1)
 
1038
    ck1 = eps / Num.sqrt(10**(0.1*rs)-1)
 
1039
    ck1p = Num.sqrt(1-ck1*ck1)
 
1040
    if ck1p == 1:
 
1041
        raise ValueError, "Cannot design a filter with given rp and rs specifications."
 
1042
 
 
1043
    wp = 1
 
1044
    val = special.ellipk([ck1*ck1,ck1p*ck1p])
 
1045
    if abs(1-ck1p*ck1p) < EPSILON:
 
1046
        krat = 0
 
1047
    else:
 
1048
        krat = N*val[0] / val[1]
 
1049
 
 
1050
    m = optimize.fmin(kratio, 0.5, args=(krat,), maxfun=250, maxiter=250,
 
1051
                      disp=0)
 
1052
    if m < 0 or m > 1:
 
1053
        m = optimize.fminbound(kratio, 0, 1, args=(krat,), maxfun=250,
 
1054
                               maxiter=250, disp=0)
 
1055
    
 
1056
    capk = special.ellipk(m)
 
1057
    ws = wp / sqrt(m)
 
1058
    m1 = 1-m
 
1059
 
 
1060
    j = Num.arange(1-N%2,N,2)
 
1061
    jj = len(j)
 
1062
 
 
1063
    [s,c,d,phi] = special.ellipj(j*capk/N,m*Num.ones(jj))
 
1064
    snew = Num.compress(abs(s) > EPSILON, s)
 
1065
    z = 1.0 / (sqrt(m)*snew)
 
1066
    z = 1j*z
 
1067
    z = Num.concatenate((z,conjugate(z)))
 
1068
 
 
1069
    r = optimize.fmin(vratio, special.ellipk(m), args=(1./eps, ck1p*ck1p),
 
1070
                      maxfun=250, maxiter=250, disp=0)
 
1071
    v0 = capk * r / (N*val[0])
 
1072
 
 
1073
    [sv,cv,dv,phi] = special.ellipj(v0,1-m)
 
1074
    p = -(c*d*sv*cv + 1j*s*dv) / (1-(d*sv)**2.0)
 
1075
 
 
1076
    if N % 2:
 
1077
        newp = Num.compress(abs(p.imag) > EPSILON*Num.sqrt(MLab.sum(p*Num.conjugate(p)).real), p)
 
1078
        p = Num.concatenate((p,conjugate(newp)))
 
1079
    else:
 
1080
        p = Num.concatenate((p,conjugate(p)))
 
1081
 
 
1082
    k = (MLab.prod(-p) / MLab.prod(-z)).real
 
1083
    if N % 2 == 0:
 
1084
        k = k / Num.sqrt((1+eps*eps))
 
1085
 
 
1086
    return z, p, k
 
1087
 
 
1088
 
 
1089
def besselap(N):
 
1090
    """Return (z,p,k) zero, pole, gain for analog prototype of an Nth
 
1091
    order Bessel filter."""
 
1092
    z = []
 
1093
    k = 1
 
1094
    if N == 0:
 
1095
        p = [];
 
1096
    elif N == 1:
 
1097
        p = [-1]
 
1098
    elif N == 2:
 
1099
        p = [-.8660254037844386467637229+.4999999999999999999999996*1j,
 
1100
             -.8660254037844386467637229-.4999999999999999999999996*1j]
 
1101
    elif N == 3:
 
1102
        p = [-.9416000265332067855971980,
 
1103
             -.7456403858480766441810907-.7113666249728352680992154*1j,
 
1104
             -.7456403858480766441810907+.7113666249728352680992154*1j]
 
1105
    elif N == 4:
 
1106
        p = [-.6572111716718829545787781-.8301614350048733772399715*1j,
 
1107
             -.6572111716718829545787788+.8301614350048733772399715*1j,
 
1108
             -.9047587967882449459642637-.2709187330038746636700923*1j,
 
1109
             -.9047587967882449459642624+.2709187330038746636700926*1j]
 
1110
    elif N == 5:
 
1111
        p = [-.9264420773877602247196260,
 
1112
             -.8515536193688395541722677-.4427174639443327209850002*1j,
 
1113
             -.8515536193688395541722677+.4427174639443327209850002*1j,
 
1114
             -.5905759446119191779319432-.9072067564574549539291747*1j,
 
1115
             -.5905759446119191779319432+.9072067564574549539291747*1j]
 
1116
    elif N == 6:
 
1117
        p = [-.9093906830472271808050953-.1856964396793046769246397*1j,
 
1118
             -.9093906830472271808050953+.1856964396793046769246397*1j,
 
1119
             -.7996541858328288520243325-.5621717346937317988594118*1j,
 
1120
             -.7996541858328288520243325+.5621717346937317988594118*1j,
 
1121
             -.5385526816693109683073792-.9616876881954277199245657*1j,
 
1122
             -.5385526816693109683073792+.9616876881954277199245657*1j]
 
1123
    elif N == 7:
 
1124
        p = [-.9194871556490290014311619,
 
1125
             -.8800029341523374639772340-.3216652762307739398381830*1j,
 
1126
             -.8800029341523374639772340+.3216652762307739398381830*1j,
 
1127
             -.7527355434093214462291616-.6504696305522550699212995*1j,
 
1128
             -.7527355434093214462291616+.6504696305522550699212995*1j,
 
1129
             -.4966917256672316755024763-1.002508508454420401230220*1j,
 
1130
             -.4966917256672316755024763+1.002508508454420401230220*1j]
 
1131
    elif N == 8:
 
1132
        p = [-.9096831546652910216327629-.1412437976671422927888150*1j,
 
1133
             -.9096831546652910216327629+.1412437976671422927888150*1j,
 
1134
             -.8473250802359334320103023-.4259017538272934994996429*1j,
 
1135
             -.8473250802359334320103023+.4259017538272934994996429*1j,
 
1136
             -.7111381808485399250796172-.7186517314108401705762571*1j,
 
1137
             -.7111381808485399250796172+.7186517314108401705762571*1j,
 
1138
             -.4621740412532122027072175-1.034388681126901058116589*1j,
 
1139
             -.4621740412532122027072175+1.034388681126901058116589*1j]
 
1140
    elif N == 9:
 
1141
        p = [-.9154957797499037686769223,
 
1142
             -.8911217017079759323183848-.2526580934582164192308115*1j,
 
1143
             -.8911217017079759323183848+.2526580934582164192308115*1j,
 
1144
             -.8148021112269012975514135-.5085815689631499483745341*1j,
 
1145
             -.8148021112269012975514135+.5085815689631499483745341*1j,
 
1146
             -.6743622686854761980403401-.7730546212691183706919682*1j,
 
1147
             -.6743622686854761980403401+.7730546212691183706919682*1j,
 
1148
             -.4331415561553618854685942-1.060073670135929666774323*1j,
 
1149
             -.4331415561553618854685942+1.060073670135929666774323*1j]
 
1150
    elif N == 10:
 
1151
        p = [-.9091347320900502436826431-.1139583137335511169927714*1j,
 
1152
             -.9091347320900502436826431+.1139583137335511169927714*1j,
 
1153
             -.8688459641284764527921864-.3430008233766309973110589*1j,
 
1154
             -.8688459641284764527921864+.3430008233766309973110589*1j,
 
1155
             -.7837694413101441082655890-.5759147538499947070009852*1j,
 
1156
             -.7837694413101441082655890+.5759147538499947070009852*1j,
 
1157
             -.6417513866988316136190854-.8175836167191017226233947*1j,
 
1158
             -.6417513866988316136190854+.8175836167191017226233947*1j,
 
1159
             -.4083220732868861566219785-1.081274842819124562037210*1j,
 
1160
             -.4083220732868861566219785+1.081274842819124562037210*1j]
 
1161
    elif N == 11:
 
1162
        p = [-.9129067244518981934637318,
 
1163
             -.8963656705721166099815744-.2080480375071031919692341*1j
 
1164
             -.8963656705721166099815744+.2080480375071031919692341*1j,
 
1165
             -.8453044014712962954184557-.4178696917801248292797448*1j,
 
1166
             -.8453044014712962954184557+.4178696917801248292797448*1j,
 
1167
             -.7546938934722303128102142-.6319150050721846494520941*1j,
 
1168
             -.7546938934722303128102142+.6319150050721846494520941*1j,
 
1169
             -.6126871554915194054182909-.8547813893314764631518509*1j,
 
1170
             -.6126871554915194054182909+.8547813893314764631518509*1j,
 
1171
             -.3868149510055090879155425-1.099117466763120928733632*1j,
 
1172
             -.3868149510055090879155425+1.099117466763120928733632*1j]
 
1173
    elif N == 12:
 
1174
        p = [-.9084478234140682638817772-95506365213450398415258360.0e-27*1j,
 
1175
             -.9084478234140682638817772+95506365213450398415258360.0e-27*1j,
 
1176
             -.8802534342016826507901575-.2871779503524226723615457*1j,
 
1177
             -.8802534342016826507901575+.2871779503524226723615457*1j,
 
1178
             -.8217296939939077285792834-.4810212115100676440620548*1j,
 
1179
             -.8217296939939077285792834+.4810212115100676440620548*1j,
 
1180
             -.7276681615395159454547013-.6792961178764694160048987*1j,
 
1181
             -.7276681615395159454547013+.6792961178764694160048987*1j,
 
1182
             -.5866369321861477207528215-.8863772751320727026622149*1j,
 
1183
             -.5866369321861477207528215+.8863772751320727026622149*1j,
 
1184
             -.3679640085526312839425808-1.114373575641546257595657*1j,
 
1185
             -.3679640085526312839425808+1.114373575641546257595657*1j]
 
1186
    elif N == 13:
 
1187
        p = [-.9110914665984182781070663,
 
1188
             -.8991314665475196220910718-.1768342956161043620980863*1j,
 
1189
             -.8991314665475196220910718+.1768342956161043620980863*1j,
 
1190
             -.8625094198260548711573628-.3547413731172988997754038*1j,
 
1191
             -.8625094198260548711573628+.3547413731172988997754038*1j,
 
1192
             -.7987460692470972510394686-.5350752120696801938272504*1j,
 
1193
             -.7987460692470972510394686+.5350752120696801938272504*1j,
 
1194
             -.7026234675721275653944062-.7199611890171304131266374*1j,
 
1195
             -.7026234675721275653944062+.7199611890171304131266374*1j,
 
1196
             -.5631559842430199266325818-.9135900338325109684927731*1j,
 
1197
             -.5631559842430199266325818+.9135900338325109684927731*1j,
 
1198
             -.3512792323389821669401925-1.127591548317705678613239*1j,
 
1199
             -.3512792323389821669401925+1.127591548317705678613239*1j]
 
1200
    elif N == 14:
 
1201
        p = [-.9077932138396487614720659-82196399419401501888968130.0e-27*1j,
 
1202
             -.9077932138396487614720659+82196399419401501888968130.0e-27*1j,
 
1203
             -.8869506674916445312089167-.2470079178765333183201435*1j,
 
1204
             -.8869506674916445312089167+.2470079178765333183201435*1j,
 
1205
             -.8441199160909851197897667-.4131653825102692595237260*1j,
 
1206
             -.8441199160909851197897667+.4131653825102692595237260*1j,
 
1207
             -.7766591387063623897344648-.5819170677377608590492434*1j,
 
1208
             -.7766591387063623897344648+.5819170677377608590492434*1j,
 
1209
             -.6794256425119233117869491-.7552857305042033418417492*1j,
 
1210
             -.6794256425119233117869491+.7552857305042033418417492*1j,
 
1211
             -.5418766775112297376541293-.9373043683516919569183099*1j,
 
1212
             -.5418766775112297376541293+.9373043683516919569183099*1j,
 
1213
             -.3363868224902037330610040-1.139172297839859991370924*1j,
 
1214
             -.3363868224902037330610040+1.139172297839859991370924*1j]
 
1215
    elif N == 15:
 
1216
        p = [-.9097482363849064167228581,
 
1217
             -.9006981694176978324932918-.1537681197278439351298882*1j,
 
1218
             -.9006981694176978324932918+.1537681197278439351298882*1j,
 
1219
             -.8731264620834984978337843-.3082352470564267657715883*1j,
 
1220
             -.8731264620834984978337843+.3082352470564267657715883*1j,
 
1221
             -.8256631452587146506294553-.4642348752734325631275134*1j,
 
1222
             -.8256631452587146506294553+.4642348752734325631275134*1j,
 
1223
             -.7556027168970728127850416-.6229396358758267198938604*1j,
 
1224
             -.7556027168970728127850416+.6229396358758267198938604*1j,
 
1225
             -.6579196593110998676999362-.7862895503722515897065645*1j,
 
1226
             -.6579196593110998676999362+.7862895503722515897065645*1j,
 
1227
             -.5224954069658330616875186-.9581787261092526478889345*1j,
 
1228
             -.5224954069658330616875186+.9581787261092526478889345*1j,
 
1229
             -.3229963059766444287113517-1.149416154583629539665297*1j,
 
1230
             -.3229963059766444287113517+1.149416154583629539665297*1j]
 
1231
    elif N == 16:
 
1232
        p = [-.9072099595087001356491337-72142113041117326028823950.0e-27*1j,
 
1233
             -.9072099595087001356491337+72142113041117326028823950.0e-27*1j,
 
1234
             -.8911723070323647674780132-.2167089659900576449410059*1j,
 
1235
             -.8911723070323647674780132+.2167089659900576449410059*1j,
 
1236
             -.8584264231521330481755780-.3621697271802065647661080*1j,
 
1237
             -.8584264231521330481755780+.3621697271802065647661080*1j,
 
1238
             -.8074790293236003885306146-.5092933751171800179676218*1j,
 
1239
             -.8074790293236003885306146+.5092933751171800179676218*1j,
 
1240
             -.7356166304713115980927279-.6591950877860393745845254*1j,
 
1241
             -.7356166304713115980927279+.6591950877860393745845254*1j,
 
1242
             -.6379502514039066715773828-.8137453537108761895522580*1j,
 
1243
             -.6379502514039066715773828+.8137453537108761895522580*1j,
 
1244
             -.5047606444424766743309967-.9767137477799090692947061*1j,
 
1245
             -.5047606444424766743309967+.9767137477799090692947061*1j,
 
1246
             -.3108782755645387813283867-1.158552841199330479412225*1j,
 
1247
             -.3108782755645387813283867+1.158552841199330479412225*1j]
 
1248
    elif N == 17:
 
1249
        p = [-.9087141161336397432860029,
 
1250
             -.9016273850787285964692844-.1360267995173024591237303*1j,
 
1251
             -.9016273850787285964692844+.1360267995173024591237303*1j,
 
1252
             -.8801100704438627158492165-.2725347156478803885651973*1j,
 
1253
             -.8801100704438627158492165+.2725347156478803885651973*1j,
 
1254
             -.8433414495836129204455491-.4100759282910021624185986*1j,
 
1255
             -.8433414495836129204455491+.4100759282910021624185986*1j,
 
1256
             -.7897644147799708220288138-.5493724405281088674296232*1j,
 
1257
             -.7897644147799708220288138+.5493724405281088674296232*1j,
 
1258
             -.7166893842372349049842743-.6914936286393609433305754*1j,
 
1259
             -.7166893842372349049842743+.6914936286393609433305754*1j,
 
1260
             -.6193710717342144521602448-.8382497252826992979368621*1j,
 
1261
             -.6193710717342144521602448+.8382497252826992979368621*1j,
 
1262
             -.4884629337672704194973683-.9932971956316781632345466*1j,
 
1263
             -.4884629337672704194973683+.9932971956316781632345466*1j,
 
1264
             -.2998489459990082015466971-1.166761272925668786676672*1j,
 
1265
             -.2998489459990082015466971+1.166761272925668786676672*1j]
 
1266
    elif N == 18:
 
1267
        p = [-.9067004324162775554189031-64279241063930693839360680.0e-27*1j,
 
1268
             -.9067004324162775554189031+64279241063930693839360680.0e-27*1j,
 
1269
             -.8939764278132455733032155-.1930374640894758606940586*1j,
 
1270
             -.8939764278132455733032155+.1930374640894758606940586*1j,
 
1271
             -.8681095503628830078317207-.3224204925163257604931634*1j,
 
1272
             -.8681095503628830078317207+.3224204925163257604931634*1j,
 
1273
             -.8281885016242836608829018-.4529385697815916950149364*1j,
 
1274
             -.8281885016242836608829018+.4529385697815916950149364*1j,
 
1275
             -.7726285030739558780127746-.5852778162086640620016316*1j,
 
1276
             -.7726285030739558780127746+.5852778162086640620016316*1j,
 
1277
             -.6987821445005273020051878-.7204696509726630531663123*1j,
 
1278
             -.6987821445005273020051878+.7204696509726630531663123*1j,
 
1279
             -.6020482668090644386627299-.8602708961893664447167418*1j,
 
1280
             -.6020482668090644386627299+.8602708961893664447167418*1j,
 
1281
             -.4734268069916151511140032-1.008234300314801077034158*1j,
 
1282
             -.4734268069916151511140032+1.008234300314801077034158*1j,
 
1283
             -.2897592029880489845789953-1.174183010600059128532230*1j,
 
1284
             -.2897592029880489845789953+1.174183010600059128532230*1j]
 
1285
    elif N == 19:
 
1286
        p = [-.9078934217899404528985092,
 
1287
             -.9021937639390660668922536-.1219568381872026517578164*1j,
 
1288
             -.9021937639390660668922536+.1219568381872026517578164*1j,
 
1289
             -.8849290585034385274001112-.2442590757549818229026280*1j,
 
1290
             -.8849290585034385274001112+.2442590757549818229026280*1j,
 
1291
             -.8555768765618421591093993-.3672925896399872304734923*1j,
 
1292
             -.8555768765618421591093993+.3672925896399872304734923*1j,
 
1293
             -.8131725551578197705476160-.4915365035562459055630005*1j,
 
1294
             -.8131725551578197705476160+.4915365035562459055630005*1j,
 
1295
             -.7561260971541629355231897-.6176483917970178919174173*1j,
 
1296
             -.7561260971541629355231897+.6176483917970178919174173*1j,
 
1297
             -.6818424412912442033411634-.7466272357947761283262338*1j,
 
1298
             -.6818424412912442033411634+.7466272357947761283262338*1j,
 
1299
             -.5858613321217832644813602-.8801817131014566284786759*1j,
 
1300
             -.5858613321217832644813602+.8801817131014566284786759*1j,
 
1301
             -.4595043449730988600785456-1.021768776912671221830298*1j,
 
1302
             -.4595043449730988600785456+1.021768776912671221830298*1j,
 
1303
             -.2804866851439370027628724-1.180931628453291873626003*1j,
 
1304
             -.2804866851439370027628724+1.180931628453291873626003*1j]
 
1305
    elif N == 20:
 
1306
        p = [-.9062570115576771146523497-57961780277849516990208850.0e-27*1j,
 
1307
             -.9062570115576771146523497+57961780277849516990208850.0e-27*1j,
 
1308
             -.8959150941925768608568248-.1740317175918705058595844*1j,
 
1309
             -.8959150941925768608568248+.1740317175918705058595844*1j,
 
1310
             -.8749560316673332850673214-.2905559296567908031706902*1j,
 
1311
             -.8749560316673332850673214+.2905559296567908031706902*1j,
 
1312
             -.8427907479956670633544106-.4078917326291934082132821*1j,
 
1313
             -.8427907479956670633544106+.4078917326291934082132821*1j,
 
1314
             -.7984251191290606875799876-.5264942388817132427317659*1j,
 
1315
             -.7984251191290606875799876+.5264942388817132427317659*1j,
 
1316
             -.7402780309646768991232610-.6469975237605228320268752*1j,
 
1317
             -.7402780309646768991232610+.6469975237605228320268752*1j,
 
1318
             -.6658120544829934193890626-.7703721701100763015154510*1j,
 
1319
             -.6658120544829934193890626+.7703721701100763015154510*1j,
 
1320
             -.5707026806915714094398061-.8982829066468255593407161*1j,
 
1321
             -.5707026806915714094398061+.8982829066468255593407161*1j,
 
1322
             -.4465700698205149555701841-1.034097702560842962315411*1j,
 
1323
             -.4465700698205149555701841+1.034097702560842962315411*1j,
 
1324
             -.2719299580251652601727704-1.187099379810885886139638*1j,
 
1325
             -.2719299580251652601727704+1.187099379810885886139638*1j]
 
1326
    elif N == 21:
 
1327
        p = [-.9072262653142957028884077,
 
1328
             -.9025428073192696303995083-.1105252572789856480992275*1j,
 
1329
             -.9025428073192696303995083+.1105252572789856480992275*1j,
 
1330
             -.8883808106664449854431605-.2213069215084350419975358*1j,
 
1331
             -.8883808106664449854431605+.2213069215084350419975358*1j,
 
1332
             -.8643915813643204553970169-.3326258512522187083009453*1j,
 
1333
             -.8643915813643204553970169+.3326258512522187083009453*1j,
 
1334
             -.8299435470674444100273463-.4448177739407956609694059*1j,
 
1335
             -.8299435470674444100273463+.4448177739407956609694059*1j,
 
1336
             -.7840287980408341576100581-.5583186348022854707564856*1j,
 
1337
             -.7840287980408341576100581+.5583186348022854707564856*1j,
 
1338
             -.7250839687106612822281339-.6737426063024382240549898*1j,
 
1339
             -.7250839687106612822281339+.6737426063024382240549898*1j,
 
1340
             -.6506315378609463397807996-.7920349342629491368548074*1j,
 
1341
             -.6506315378609463397807996+.7920349342629491368548074*1j,
 
1342
             -.5564766488918562465935297-.9148198405846724121600860*1j,
 
1343
             -.5564766488918562465935297+.9148198405846724121600860*1j,
 
1344
             -.4345168906815271799687308-1.045382255856986531461592*1j,
 
1345
             -.4345168906815271799687308+1.045382255856986531461592*1j,
 
1346
             -.2640041595834031147954813-1.192762031948052470183960*1j,
 
1347
             -.2640041595834031147954813+1.192762031948052470183960*1j]
 
1348
    elif N == 22:
 
1349
        p = [-.9058702269930872551848625-52774908289999045189007100.0e-27*1j,
 
1350
             -.9058702269930872551848625+52774908289999045189007100.0e-27*1j,
 
1351
             -.8972983138153530955952835-.1584351912289865608659759*1j,
 
1352
             -.8972983138153530955952835+.1584351912289865608659759*1j,
 
1353
             -.8799661455640176154025352-.2644363039201535049656450*1j,
 
1354
             -.8799661455640176154025352+.2644363039201535049656450*1j,
 
1355
             -.8534754036851687233084587-.3710389319482319823405321*1j,
 
1356
             -.8534754036851687233084587+.3710389319482319823405321*1j,
 
1357
             -.8171682088462720394344996-.4785619492202780899653575*1j,
 
1358
             -.8171682088462720394344996+.4785619492202780899653575*1j,
 
1359
             -.7700332930556816872932937-.5874255426351153211965601*1j,
 
1360
             -.7700332930556816872932937+.5874255426351153211965601*1j,
 
1361
             -.7105305456418785989070935-.6982266265924524000098548*1j,
 
1362
             -.7105305456418785989070935+.6982266265924524000098548*1j,
 
1363
             -.6362427683267827226840153-.8118875040246347267248508*1j,
 
1364
             -.6362427683267827226840153+.8118875040246347267248508*1j,
 
1365
             -.5430983056306302779658129-.9299947824439872998916657*1j,
 
1366
             -.5430983056306302779658129+.9299947824439872998916657*1j,
 
1367
             -.4232528745642628461715044-1.055755605227545931204656*1j,
 
1368
             -.4232528745642628461715044+1.055755605227545931204656*1j,
 
1369
             -.2566376987939318038016012-1.197982433555213008346532*1j,
 
1370
             -.2566376987939318038016012+1.197982433555213008346532*1j]
 
1371
    elif N == 23:
 
1372
        p = [-.9066732476324988168207439,
 
1373
             -.9027564979912504609412993-.1010534335314045013252480*1j,
 
1374
             -.9027564979912504609412993+.1010534335314045013252480*1j,
 
1375
             -.8909283242471251458653994-.2023024699381223418195228*1j,
 
1376
             -.8909283242471251458653994+.2023024699381223418195228*1j,
 
1377
             -.8709469395587416239596874-.3039581993950041588888925*1j,
 
1378
             -.8709469395587416239596874+.3039581993950041588888925*1j,
 
1379
             -.8423805948021127057054288-.4062657948237602726779246*1j,
 
1380
             -.8423805948021127057054288+.4062657948237602726779246*1j,
 
1381
             -.8045561642053176205623187-.5095305912227258268309528*1j,
 
1382
             -.8045561642053176205623187+.5095305912227258268309528*1j,
 
1383
             -.7564660146829880581478138-.6141594859476032127216463*1j,
 
1384
             -.7564660146829880581478138+.6141594859476032127216463*1j,
 
1385
             -.6965966033912705387505040-.7207341374753046970247055*1j,
 
1386
             -.6965966033912705387505040+.7207341374753046970247055*1j,
 
1387
             -.6225903228771341778273152-.8301558302812980678845563*1j,
 
1388
             -.6225903228771341778273152+.8301558302812980678845563*1j,
 
1389
             -.5304922463810191698502226-.9439760364018300083750242*1j,
 
1390
             -.5304922463810191698502226+.9439760364018300083750242*1j,
 
1391
             -.4126986617510148836149955-1.065328794475513585531053*1j,
 
1392
             -.4126986617510148836149955+1.065328794475513585531053*1j,
 
1393
             -.2497697202208956030229911-1.202813187870697831365338*1j,
 
1394
             -.2497697202208956030229911+1.202813187870697831365338*1j]
 
1395
    elif N == 24:
 
1396
        p = [-.9055312363372773709269407-48440066540478700874836350.0e-27*1j,
 
1397
             -.9055312363372773709269407+48440066540478700874836350.0e-27*1j,
 
1398
             -.8983105104397872954053307-.1454056133873610120105857*1j,
 
1399
             -.8983105104397872954053307+.1454056133873610120105857*1j,
 
1400
             -.8837358034555706623131950-.2426335234401383076544239*1j,
 
1401
             -.8837358034555706623131950+.2426335234401383076544239*1j,
 
1402
             -.8615278304016353651120610-.3403202112618624773397257*1j,
 
1403
             -.8615278304016353651120610+.3403202112618624773397257*1j,
 
1404
             -.8312326466813240652679563-.4386985933597305434577492*1j,
 
1405
             -.8312326466813240652679563+.4386985933597305434577492*1j,
 
1406
             -.7921695462343492518845446-.5380628490968016700338001*1j,
 
1407
             -.7921695462343492518845446+.5380628490968016700338001*1j,
 
1408
             -.7433392285088529449175873-.6388084216222567930378296*1j,
 
1409
             -.7433392285088529449175873+.6388084216222567930378296*1j,
 
1410
             -.6832565803536521302816011-.7415032695091650806797753*1j,
 
1411
             -.6832565803536521302816011+.7415032695091650806797753*1j,
 
1412
             -.6096221567378335562589532-.8470292433077202380020454*1j,
 
1413
             -.6096221567378335562589532+.8470292433077202380020454*1j,
 
1414
             -.5185914574820317343536707-.9569048385259054576937721*1j,
 
1415
             -.5185914574820317343536707+.9569048385259054576937721*1j,
 
1416
             -.4027853855197518014786978-1.074195196518674765143729*1j,
 
1417
             -.4027853855197518014786978+1.074195196518674765143729*1j,
 
1418
             -.2433481337524869675825448-1.207298683731972524975429*1j,
 
1419
             -.2433481337524869675825448+1.207298683731972524975429*1j]
 
1420
    elif N == 25:
 
1421
        p = [-.9062073871811708652496104,
 
1422
             -.9028833390228020537142561-93077131185102967450643820.0e-27*1j,
 
1423
             -.9028833390228020537142561+93077131185102967450643820.0e-27*1j,
 
1424
             -.8928551459883548836774529-.1863068969804300712287138*1j,
 
1425
             -.8928551459883548836774529+.1863068969804300712287138*1j,
 
1426
             -.8759497989677857803656239-.2798521321771408719327250*1j,
 
1427
             -.8759497989677857803656239+.2798521321771408719327250*1j,
 
1428
             -.8518616886554019782346493-.3738977875907595009446142*1j,
 
1429
             -.8518616886554019782346493+.3738977875907595009446142*1j,
 
1430
             -.8201226043936880253962552-.4686668574656966589020580*1j,
 
1431
             -.8201226043936880253962552+.4686668574656966589020580*1j,
 
1432
             -.7800496278186497225905443-.5644441210349710332887354*1j,
 
1433
             -.7800496278186497225905443+.5644441210349710332887354*1j,
 
1434
             -.7306549271849967721596735-.6616149647357748681460822*1j,
 
1435
             -.7306549271849967721596735+.6616149647357748681460822*1j,
 
1436
             -.6704827128029559528610523-.7607348858167839877987008*1j,
 
1437
             -.6704827128029559528610523+.7607348858167839877987008*1j,
 
1438
             -.5972898661335557242320528-.8626676330388028512598538*1j,
 
1439
             -.5972898661335557242320528+.8626676330388028512598538*1j,
 
1440
             -.5073362861078468845461362-.9689006305344868494672405*1j,
 
1441
             -.5073362861078468845461362+.9689006305344868494672405*1j,
 
1442
             -.3934529878191079606023847-1.082433927173831581956863*1j,
 
1443
             -.3934529878191079606023847+1.082433927173831581956863*1j,
 
1444
             -.2373280669322028974199184-1.211476658382565356579418*1j,
 
1445
             -.2373280669322028974199184+1.211476658382565356579418*1j]
 
1446
    else:
 
1447
        raise ValueError, "Bessel Filter not supported for order %d" % N
 
1448
 
 
1449
    return z, p, k
 
1450
 
 
1451
filter_dict = {'butter': [buttap,buttord],
 
1452
               'butterworth' : [buttap,buttord],
 
1453
               'cauer' : [ellipap,ellipord],
 
1454
               'elliptic' : [ellipap,ellipord],
 
1455
               'ellip' : [ellipap,ellipord],
 
1456
               'bessel' : [besselap],
 
1457
               'cheby1' : [cheb1ap, cheb1ord],
 
1458
               'chebyshev1' : [cheb1ap, cheb1ord],
 
1459
               'chebyshevi' : [cheb1ap, cheb1ord],
 
1460
               'cheby2' : [cheb2ap, cheb2ord],
 
1461
               'chebyshev2' : [cheb2ap, cheb2ord],
 
1462
               'chebyshevii' : [cheb2ap, cheb2ord]
 
1463
               }
 
1464
 
 
1465
band_dict = {'band':'bandpass',
 
1466
             'bandpass':'bandpass',
 
1467
             'pass' : 'bandpass',
 
1468
             'bp':'bandpass',
 
1469
             'bs':'bandstop',
 
1470
             'bandstop':'bandstop',
 
1471
             'bands' : 'bandstop',
 
1472
             'stop' : 'bandstop',
 
1473
             'l' : 'lowpass',
 
1474
             'low': 'lowpass',
 
1475
             'lowpass' : 'lowpass',
 
1476
             'high' : 'highpass',
 
1477
             'highpass' : 'highpass',
 
1478
             'h' : 'highpass'
 
1479
             }
 
1480
             
 
1481
def kaiserord(ripple, width):
 
1482
    """Design a Kaiser window to limit ripple and width of transition region.
 
1483
 
 
1484
    Inputs:
 
1485
 
 
1486
      ripple -- positive number specifying maximum ripple in passband (dB)
 
1487
                  and minimum ripple in stopband
 
1488
      width  -- width of transition region (normalized so that 1 corresponds
 
1489
                  to pi radians / sample)
 
1490
 
 
1491
    Outputs:
 
1492
 
 
1493
      N, beta -- the order and beta parameter for the kaiser window.
 
1494
 
 
1495
                 signal.kaiser(N,beta,sym=0) returns the window as does
 
1496
                 signal.get_window(beta,N)
 
1497
                 signal.get_window(('kaiser',beta),N)
 
1498
 
 
1499
    Uses the empirical equations discovered by Kaiser.
 
1500
 
 
1501
    Oppenheim, Schafer, "Discrete-Time Signal Processing,", p.475-476.
 
1502
    """
 
1503
    A = abs(ripple)  # in case somebody is confused as to what's meant
 
1504
    beta = select([A>50, A>21],
 
1505
                  [0.1102*(A-8.7), 0.5842*(A-21)**(0.4) + 0.07866*(A-21)],
 
1506
                  0.0)
 
1507
    N = (A-8)/2.285/(pi*width)
 
1508
    return ceil(N), beta 
 
1509
 
 
1510
def firwin(N, cutoff, width=None, window='hamming'):
 
1511
    """FIR Filter Design using windowed ideal filter method.
 
1512
 
 
1513
    Inputs:
 
1514
    
 
1515
      N      -- order of filter (number of taps)
 
1516
      cutoff -- cutoff frequency of filter (normalized so that 1 corresponds to
 
1517
                  Nyquist or pi radians / sample)
 
1518
 
 
1519
      width  -- if width is not None, then assume it is the approximate width of
 
1520
                  the transition region (normalized so that 1 corresonds to pi)
 
1521
                  for use in kaiser FIR filter design.
 
1522
      window -- desired window to use.
 
1523
 
 
1524
    Outputs:
 
1525
 
 
1526
      h      -- coefficients of length N fir filter. 
 
1527
    """
 
1528
 
 
1529
    from signaltools import get_window
 
1530
    if isinstance(width,float):
 
1531
        A = 2.285*N*width + 8
 
1532
        if (A < 21): beta = 0.0
 
1533
        elif (A <= 50): beta = 0.5842*(A-21)**0.4 + 0.07886*(A-21)
 
1534
        else: beta = 0.1102*(A-8.7)
 
1535
        window=('kaiser',beta)
 
1536
 
 
1537
    win = get_window(window,N,fftbins=1)
 
1538
    alpha = N//2
 
1539
    m = Num.arange(0,N)
 
1540
    h = win*special.sinc(cutoff*(m-alpha))
 
1541
    return h / sum(h)
 
1542
    
 
1543
    
 
1544