```python
from scipy import signal
from bokeh.plotting import figure, output_file, show
from bokeh.io import output_notebook
from bokeh.layouts import row
import numpy as np
import math
import sympy as sy
output_notebook()
```



<div class="bk-root">
    <a href="http://bokeh.pydata.org" target="_blank" class="bk-logo bk-logo-small bk-logo-notebook"></a>
    <span id="6e81394f-6e1b-4d92-9353-24799e182188">Loading BokehJS ...</span>
</div>





```python
def lti_to_sympy(lsys, symplify=True):
    """ Convert Scipy's LTI instance to Sympy expression """
    s = sy.Symbol('s')
    G = sy.Poly(lsys.num, s) / sy.Poly(lsys.den, s)
    return sy.simplify(G) if symplify else G


def sympy_to_lti(xpr, s=sy.Symbol('s')):
    """ Convert Sympy transfer function polynomial to Scipy LTI """
    num, den = sy.simplify(xpr).as_numer_denom()  # expressions
    p_num_den = sy.poly(num, s), sy.poly(den, s)  # polynomials
    c_num_den = [sy.expand(p).all_coeffs() for p in p_num_den]  # coefficients
    l_num, l_den = [sy.lambdify((), c)() for c in c_num_den]  # convert to floats
    return signal.lti(l_num, l_den)


def polyphase_iir_to_ba_coeffs(directCoeffs,delayedCoeffs):
    #For z^-1:
    z_1 = sy.Symbol('z_1')
    directPath = 1
    delayedPath = 1
    for c in directCoeffs:
        directPath = directPath * (sy.Poly([1,0,c], z_1) / sy.Poly([c,0,1], z_1))
    for c in delayedCoeffs:
        delayedPath = delayedPath * (sy.Poly([1,0,c], z_1) / sy.Poly([c,0,1], z_1))
        
    tranferFunction = 0.5 *( directPath + sy.Poly([1,0], z_1) * delayedPath)
    n,d = sy.fraction(sy.simplify(tranferFunction))
    b = [float(c) for c in sy.Poly(n).all_coeffs()]
    a = [float(c) for c in sy.Poly(d).all_coeffs()]
    return b,a

""" Adapted from:
    Digital signal processing schemes for efficient interpolation and decimation
    R.A. Velenzuela, A.G. Constantinides
    IEEE 1983
    
    transitionBandwitdth between 0 and 0.5 ( 0.5 = half nyquist)
    stopbandAttenuation in dB
"""
def design_polyphase_halfband_iir(transitionBandwidth, stopbandAttenuation):
    k = np.tan((np.pi-2*np.pi*transitionBandwidth)/4)**2
    kp = np.sqrt(1-k**2)
    e = 0.5 * (1 -np.sqrt(kp))/(1+np.sqrt(kp))
    q = e + 2*e**5 + 15*e**9 + 150*e**13
    ds = 10**(-stopbandAttenuation/20)
    k1 = ds**2 / (1 - ds**2)
    n = int(math.ceil(math.log(k1**2/16)/math.log(q)))
    if n % 2 ==0:
        n += 1
    if n ==1 :
        n = 3
    print("Order: %d" % n)
    q1 = q**n
    k1 = 4 * math.sqrt(q1)
    ds = math.sqrt(k1 / (1+k1))
    dp = 1 - math.sqrt(1-ds**2)
    
    def ellipticSum1(q,n,i):
        s = 0
        for m in range(5):
            s += (-1)**m *q**(m*(m+1)) * math.sin((2*m+1)*(math.pi*i)/n)
        return s
    def ellipticSum2(q,n,i):
        s = 0
        for m in range(1,5):
            s += (-1)**m *q**(m*m) * math.cos(2*m*(math.pi*i)/n)
        return s
    wi = [ 2*math.pow(q,0.25) * ellipticSum1(q,n,i)/(1+2*ellipticSum2(q,n,i)) for i in range(1,int((n-1)/2)+1) ]
    ai = [ math.sqrt((1-(w**2)*k)*(1-(w**2)/k))/(1+w**2) for w in wi ]
    ai = [ (1-a)/(1+a) for a in ai ]
    #print(ai)
    return ai[0::2], ai[1::2]
        
```


```python
def plot_filter(b,a,name):
    p1 = figure(plot_width = 500, plot_height=300, title = "%s Gain" % name, y_range = (-80,6))
    p1.xaxis.axis_label = 'Frequency [rad/sample]'
    p1.yaxis.axis_label = 'Amplitude [dB]'
    w, h = signal.freqz(b,a)
    p1.line(w, 20 * np.log10(abs(h)))
    
    w, gd = signal.group_delay((b, a))
    p2 = figure(plot_width = 500, plot_height=300, title = "%s Group delay" % name)
    p2.xaxis.axis_label = 'Frequency [rad/sample]'
    p2.yaxis.axis_label = 'Group delay [samples]'
    p2.line(w,gd)
    show(row(p1,p2))
```


```python
#b, a = signal.iirdesign(0., 0.3, 5, 50, ftype='cheby1')
b,a  = signal.cheby1(5,1,0.5)
np.set_printoptions(precision=16)
#print(a)
#print(b)
plot_filter(b,a,'Cheby1')
```




    <div class="bk-root">
        <div class="bk-plotdiv" id="d45044c6-afc2-420f-aa81-30bc61023c93"></div>
    </div>




```python
b,a  = signal.butter(7,0.5)
plot_filter(b,a,'Butterworth')
```

    C:\Users\jmarin\AppData\Local\Continuum\Anaconda3\lib\site-packages\scipy\signal\filter_design.py:536: UserWarning: The group delay is singular at frequencies [3.129, 3.135], setting to 0
      format(", ".join("{0:.3f}".format(ws) for ws in w[singular]))





    <div class="bk-root">
        <div class="bk-plotdiv" id="3d9597cb-f14f-4b23-b61f-662646c89829"></div>
    </div>




```python
b,a  = signal.cheby2(5,40,0.5)
plot_filter(b,a,'Cheby2')
```




    <div class="bk-root">
        <div class="bk-plotdiv" id="06bfc4c9-52da-4e6f-80b5-2a992aa15b1c"></div>
    </div>




```python
b,a = polyphase_iir_to_ba_coeffs([1.0 / (5 + 2 * np.sqrt(5))],[5-2*np.sqrt(5)])
plot_filter(b,a,'ButterPolyphase')
```




    <div class="bk-root">
        <div class="bk-plotdiv" id="5f52cb6a-6331-4795-b4a1-68e646923885"></div>
    </div>




```python
int(math.ceil(0.9))
```




    1




```python
directCoeffs,delayedCoeffs = design_polyphase_halfband_iir(0.1, 30)
print(directCoeffs)
print(delayedCoeffs)
b,a = polyphase_iir_to_ba_coeffs(directCoeffs,delayedCoeffs)
plot_filter(b,a,'Cheby-ish polyphase')
```

    Order: 5
    [0.23647102099689202]
    [0.71454214971259888]





    <div class="bk-root">
        <div class="bk-plotdiv" id="06d3c9f5-bb0c-4512-a26d-8fad771d30e6"></div>
    </div>




```python
[i for i in range(1,int((7-1)/2)+1)]
```




    [1, 2, 3]




```python
[1.0 / (5 + 2 * np.sqrt(5))],[5-2*np.sqrt(5)]
```




    ([0.10557280900008412], [0.52786404500042039])




```python

```

b6a6ed0f382dedd87b9a1b7e40c289b2f6e0b319
ipynb
Jupyter Notebook
tests/FiltersTests.ipynb
chuckhousley/TriggerFish-VCV
e16db7397695d61dc067f5a31140ecfd201253eb
BSD-3-Clause
2018-09-17T01:16:29.000Z
2021-06-09T16:20:21.000Z
tests/FiltersTests.ipynb
chuckhousley/TriggerFish-VCV
e16db7397695d61dc067f5a31140ecfd201253eb
BSD-3-Clause
2018-09-16T12:51:30.000Z
2020-02-13T20:31:27.000Z
tests/FiltersTests.ipynb
chuckhousley/TriggerFish-VCV
e16db7397695d61dc067f5a31140ecfd201253eb
BSD-3-Clause
2019-08-21T20:56:41.000Z
2022-01-08T14:13:44.000Z
Qwen/Qwen-72B
1. YES
2. YES