%matplotlib inline

import numpy as np
import scipy
import scipy.signal
import matplotlib.pyplot as plt
import ipywidgets as widgets

import sympy
from sympy.interactive import printing
printing.init_printing(use_latex=True)

Analog version

Q, omega_0, s = sympy.symbols("Q omega_0 s")
H_a = -omega_0**2/(s**2 + s*(omega_0/Q) + omega_0**2)
H_a

$$- \frac{\omega_{0}^{2}}{\omega_{0}^{2} + s^{2} + \frac{\omega_{0} s}{Q}}$$

Compare frequency response

Compare the frequency response for

• analog
• rbj cookbook
• Hal chamberlin version
• My derived version

Q, omega, T, z = sympy.symbols("Q omega T z")

factor = Q*sympy.tan(T*omega/2)**2 + Q + sympy.tan(T*omega/2)
a_0 = 1
a_1 = ((2*Q*sympy.tan(T*omega/2)**2 - 2*Q)/factor)
a_2 = ((Q*sympy.tan(T*omega/2)**2 + Q - sympy.tan(T*omega/2))/factor)
b_0 = (Q*sympy.tan(T*omega/2)**2/factor)
b_1 = (2*b_0)
b_2 = (b_0)

H = (b_0*z**2 + b_1*z**1 + b_2*z**0)/(a_0*z**2 + a_1*z**1 + a_2*z**0)
H

$$\frac{\frac{Q z^{2} \tan^{2}{\left (\frac{T \omega}{2} \right )}}{Q \tan^{2}{\left (\frac{T \omega}{2} \right )} + Q + \tan{\left (\frac{T \omega}{2} \right )}} + \frac{2 Q z \tan^{2}{\left (\frac{T \omega}{2} \right )}}{Q \tan^{2}{\left (\frac{T \omega}{2} \right )} + Q + \tan{\left (\frac{T \omega}{2} \right )}} + \frac{Q \tan^{2}{\left (\frac{T \omega}{2} \right )}}{Q \tan^{2}{\left (\frac{T \omega}{2} \right )} + Q + \tan{\left (\frac{T \omega}{2} \right )}}}{z^{2} + \frac{z \left(2 Q \tan^{2}{\left (\frac{T \omega}{2} \right )} - 2 Q\right)}{Q \tan^{2}{\left (\frac{T \omega}{2} \right )} + Q + \tan{\left (\frac{T \omega}{2} \right )}} + \frac{Q \tan^{2}{\left (\frac{T \omega}{2} \right )} + Q - \tan{\left (\frac{T \omega}{2} \right )}}{Q \tan^{2}{\left (\frac{T \omega}{2} \right )} + Q + \tan{\left (\frac{T \omega}{2} \right )}}}$$

H_tmp = H

# Set wanted values
sample_freq = 44.1E3
T_val = 1/sample_freq
Q_val = 1/np.sqrt(2)

w = np.array([0])
h = np.array([0])
def f_freqz(f_cutoff,Q_val):
    omega_cutoff = f_cutoff * 2 * np.pi #in rads/sec
    omega_val = omega_cutoff
    
    #Plot analog
    a_a = np.array([1, omega_val/Q_val, omega_val**2])
    b_a = np.array([0, 0, omega_val**2])
    w_a, h_a = scipy.signal.freqs(b_a, a_a, 1000)
    w_a_hz = w_a/(2*np.pi)

    plt.subplot(4, 1, 1)
    plt.plot(w_a_hz, 20 * np.log10(abs(h_a)), 'b')
    plt.title('Frequency response analog (f0=%g Hz, Q=%g)' % (f_cutoff, Q_val), color='b')
    plt.ylabel('Amplitude [dB]', color='b')
    plt.xlabel('Frequency [Hz]')
    plt.xlim([0, sample_freq/2])
    plt.ylim([-20, 20])
    plt.show()
    
    #plot rbj
    w0 = 2*np.pi*f_cutoff/sample_freq
    alpha = np.sin(w0)/(2*Q_val)
    b_0_rbj_val = (1-np.cos(w0))/2
    b_1_rbj_val = (1-np.cos(w0))
    b_2_rbj_val = (1-np.cos(w0))/2
    a_0_rbj_val = (1+alpha)
    a_1_rbj_val = -2*np.cos(w0)
    a_2_rbj_val = 1-alpha
    
    b_rbj = np.array([b_0_rbj_val, b_1_rbj_val, b_2_rbj_val],np.float)
    a_rbj = np.array([a_0_rbj_val, a_1_rbj_val, a_2_rbj_val],np.float)
    
    w_rbj, h_rbj = scipy.signal.freqz(b_rbj, a_rbj, 100)
    plt.subplot(4, 1, 2)
    plt.plot(w_rbj, 20 * np.log10(abs(h_rbj)), 'b')
    plt.title('Frequency response digital  RBJ IIR (f0=%g Hz, Q=%g)' % (f_cutoff, Q_val), color='b')
    plt.ylabel('Amplitude [dB]', color='b')
    plt.xlabel('Normalized Frequency [rad/sample]')
    plt.ylim([-20, 20])
    plt.xlim([0, np.pi])
    plt.show()
    
    #plot HalCham
    w0 = 2*np.pi*f_cutoff/sample_freq
    kf = 2*np.sin(w0/2)
    kq = 1/Q_val
    b_0_hal_val = kf*kf
    b_1_hal_val = 0
    b_2_hal_val = 0
    a_0_hal_val = 1
    a_1_hal_val = -2+kf*(kf+kq)
    a_2_hal_val = 1 - (kf*kq)
    
    b_hal = np.array([b_0_hal_val],np.float)
    a_hal = np.array([a_0_hal_val, a_1_hal_val, a_2_hal_val],np.float)
    
    w_hal, h_hal = scipy.signal.freqz(b_hal, a_hal)
    plt.subplot(4, 1, 3)
    plt.plot(w_hal, 20 * np.log10(abs(h_hal)), 'b')
    plt.title('Frequency response digital  HALCHAM IIR (f0=%g Hz, Q=%g)' % (f_cutoff, Q_val), color='b')
    plt.ylabel('Amplitude [dB]', color='b')
    plt.xlabel('Normalized Frequency [rad/sample]')
    plt.ylim([-20, 20])
    plt.xlim([0, np.pi])
    plt.show()
    
    #plot digital
    omega_cutoff = f_cutoff * 2 * np.pi #in rads/sec
    omega_val = omega_cutoff

    b_0_val = b_0.subs(T,T_val).subs(Q,Q_val).subs(omega,omega_val)
    b_1_val = b_1.subs(T,T_val).subs(Q,Q_val).subs(omega,omega_val)
    b_2_val = b_2.subs(T,T_val).subs(Q,Q_val).subs(omega,omega_val)
    a_0_val = a_0
    a_1_val = a_1.subs(T,T_val).subs(Q,Q_val).subs(omega,omega_val)
    a_2_val = a_2.subs(T,T_val).subs(Q,Q_val).subs(omega,omega_val)
    
    b = np.array([b_0_val, b_1_val, b_2_val],np.float)
    a = np.array([a_0_val, a_1_val, a_2_val],np.float)
    
    w, h = scipy.signal.freqz(b, a, 100)
    plt.subplot(4, 1, 4)
    plt.plot(w, 20 * np.log10(abs(h)), 'b')
    plt.title('Frequency response digital IIR (f0=%g Hz, Q=%g)' % (f_cutoff, Q_val), color='b')
    plt.ylabel('Amplitude [dB]', color='b')
    plt.xlabel('Normalized Frequency [rad/sample]')
    plt.ylim([-20, 20])
    plt.xlim([0, np.pi])
    plt.show()

widgets.interact(f_freqz,f_cutoff=widgets.IntSlider(min=20,max=22.05E3,step=1,value=15E3),Q_val=widgets.FloatSlider(min=0.1,max=8,step=0.05,value=1.0))