Laplace transfer tunction for filter equivalent to Matlab 2nd order butterworth filter

Question

I am trying to build a bandpass filter analytically with the following specification:

1. Center frequency of fc (say $16Hz$)
2. Pass bandwidth of fb (say 2hz. therefore, cut-off frequencies are $15$ and $17Hz$)
3. $0\ \rm dB$ in the pass bandwidth

The goal is to develop a filter that matches the one generated by the following code:

freq_Plox = 16;% [Hz]
freq_sepa = 1; % separation width
ord       = 2; % butterworth filter order
wc_bpf_lp = 2*pi*(freq_Plox - freq_sepa);
wc_bpf_hp = 2*pi*(freq_Plox + freq_sepa);
dt = 1/1000;

[num_bpf,den_bpf] = butter(ord,(2*dt)/(2*pi)*[wc_bpf_lp wc_bpf_hp],'bandpass');
opts = bodeoptions;
opts.PhaseWrapping = 'on';
opts.FreqUnits = 'Hz';

figure;bode(tf(num_bpf,den_bpf,dt_mec_di),opts);grid on;hold on;

Without a good textbook on filters, i resorted to web search and none of the formulas has worked. Please, help.

[UPDATE]

I made some progress. The following code gets really close to the Matlab filter except that the gain is off.

fc = 16; %[Hz], center frequency
fb = 1;  %[Hz], ±pass bandwidth
fl = fc - fb; %[Hz], low cutoff frequency
fh = fc + fb; %[Hz], high cutoff frequency
Q  = fc / (fh - fl); % quality factor

% cascade 4th order
wc = 2*pi*fc; %[rad/sec]
wl = 2*pi*fl; %[rad/sec]
wh = 2*pi*fh; %[rad/sec]

Hs = Kpb* wh^2/(s^2 + (wh/Q)*s + wh^2) * s^2/(s^2 + (wl/Q)*s + wl^2);

Then, the final question is how to analytically calculate the Kpb that will set the peak amplitude to $0\ \rm dB$. When I plotted the low-pass and the high-pass filter separately, at their natural frequency, $w_l$ and $w_h$, the peak amplitudes are $\sim 18\ \rm dB$ for both; totaling $36\rm dB$. What is the formula that will result in $36\ \rm dB$?

Answer 1

Your filter transfer function in [UPDATE]:

$$H= \left(K_{\mathit{pb}},s\right)\rightarrow K_{\mathit{pb}}\cdot \frac{\omega_{h}^{2}}{s^{2}+\frac{\omega_{h}}{Q}\cdot s+\omega_{h}^{2}}\cdot \frac{s^{2}}{s^{2}+\frac{\omega_{l}}{Q}\cdot s+\omega_{l}^{2}}$$

With $s= i\cdot \omega$ we get the absolute value:

$$H_{\mathrm{abs}}= \left(K_{\mathit{pb}},\omega \right)\mapsto \frac{\omega_{h}^{2}\cdot \omega^{2}\cdot K_{\mathit{pb}}\cdot Q^{2}}{\sqrt{\left(\omega -\omega_{l}\right)^{2}\cdot \left(\omega +\omega_{l}\right)^{2}\cdot Q^{2}+\omega_{l}^{2}\cdot \omega^{2}}\cdot \sqrt{\left(\omega -\omega_{h}\right)^{2}\cdot \left(\omega +\omega_{h}\right)^{2}\cdot Q^{2}+\omega_{h}^{2}\cdot \omega^{2}}}$$

Finally we determine $K_{pb}$ so that $H_{abs}(K_{pb},\omega_c)=1$:

$$K_{\mathit{pb}}= \frac{1}{H_{\mathrm{abs}}\left(1,\omega_{c}\right)}=\frac{\sqrt{\left(\omega_{c}-\omega_{l}\right)^{2} \left(\omega_{c}+\omega_{l}\right)^{2} Q^{2}+\omega_{l}^{2} \omega_{c}^{2}}\, \sqrt{\left(\omega_{c}-\omega_{h}\right)^{2} \left(\omega_{c}+\omega_{h}\right)^{2} Q^{2}+\omega_{h}^{2} \omega_{c}^{2}}}{\omega_{h}^{2} \omega_{c}^{2} Q^{2}}$$

Alternatively because the peak amplitude is at $\omega_{peak}=\sqrt{\omega_{l} \omega_{h}}$ (geometric mean):

$$K_{\mathit{pb}}=\frac{1}{H_{\mathrm{abs}}\left(1,\sqrt{\omega_{l}\omega_{h}}\right)}= \frac{Q^{2} \omega_{h}^{2}+\omega_{l} \left(-2 Q^{2}+1\right) \omega_{h}+Q^{2} \omega_{l}^{2}}{Q^{2} \omega_{h}^{2}}$$