import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt
# Exploring a so-called Linkwitz Transform filter for subwoofer response tailoring
fs = 44100
fn = fs//2 # Nyquist
fc = 36.95 # corner freq of closed-box system in Hz
Q = 0.803 # Q of close-box system (aka Q_tc)
b = [1,0,0] # zeros
a = [1, (2*np.pi*fc/Q), (2*np.pi*fc)**2] # poles
w = np.linspace(5, 2*np.pi*fn, num=fn+1, endpoint=True)
w, Ho = sig.freqs(b, a, worN=w)
plt.figure()
plt.semilogx(w/(2*np.pi), 20*np.log10(np.abs(Ho)), label='Original System')
#plt.xlim([1, 150])
fc1 = 20 # desired corner frequency of close-box system in Hz
Q1 = 1/np.sqrt(2) # desired Q of closed-box system # np.euler_gamma #
lt_b = a
lt_a = [1, (2*np.pi*fc1/Q1), (2*np.pi*fc1)**2]
w1, H1 = sig.freqs(lt_b, lt_a, worN=w)
plt.semilogx(w1/(2*np.pi), 20*np.log10(np.abs(H1)), label='LT Filter')
bt = np.convolve(b, lt_b)
at = np.convolve(a, lt_a)
wt, Ht = sig.freqs(bt, at, worN=w)
plt.semilogx(wt/(2*np.pi), 20*np.log10(np.abs(Ht)), label='LT + Original')
plt.axhline(-3, color='grey', linestyle='--')
plt.axvline(20, color='grey', linestyle='--')
plt.xlabel('Freq [Hz]')
plt.legend(loc='lower right')
plt.ylim([-30, np.max(20*np.log10(np.abs(H1)))+2])
plt.xlim([5, 1200])
plt.grid(b=True, which='both', axis='both')