Cascade of First Order High-shelf and Low-shelf Sections
This implementation is simply the cascade of two first order sections with real coefficients and thus is limited to what can be accomplished with poles and zeros on the real axis as first-order high-shelf and low-shelf sections. It appears this was done by mapping the continuous time transfer function for simple first order high-shelf and low-shelf filters to discrete time using the Bilinear Transform along with prewarping the cutoff frequencies to compensate for the frequency warping that occurs with that mapping technique. The cascade of two first order sections in series can be combined as a single biquad filter as the product of their individual transfer functions. Here we see the resulting first-order filters as:
$$H_L(z) = \frac{\beta_1z+\beta_2}{\alpha_0z+\alpha_1}$$
$$H_H(z) = \frac{\beta_3z+\beta_4}{\alpha_2z+\alpha_3}$$
The cascade would be the product as follows:
$$H_L(z)H_H(z) = \frac{(\beta_1z+\beta_2)(\beta_3z+\beta_4)}{(\alpha_0z+\alpha_1)(\alpha_2z+\alpha_3)}$$
$$= \frac{(\beta_1\beta_3)z^2+(\beta_1 \beta_4 + \beta_2 \beta_3)z + \beta_2 \beta_4}{(\alpha_0\alpha_2) z^2+ (\alpha_0 \alpha_3 + \alpha_1 \alpha_2)z + \alpha_1\alpha_3}$$
And thus the biquad coefficients are given as:
$b_0 = \beta_1\beta_3$
$b_1 = \beta_1 \beta_4 + \beta_2 \beta_3$
$b_2 = \beta_2 \beta_4$
$a_0 = \alpha_0\alpha_2$
$a_1 = \alpha_0 \alpha_3 + \alpha_1 \alpha_2$
$a_2 = \alpha_1\alpha_3$
For example, the low-shelf (boosting the Bass) is accomplished by placing the pole closer to the origin than the zero on the real axis in the S Plane, and the graphic below illustrates how this would map to a pole and zero placement in the z-plane using the Bilinear Transform resulting in a pole at $z=-\alpha_1/\alpha_0$ and a zero at $z=-\beta_2/\beta_1$ where the origin in the s-plane maps to $z=1$, and the zero is now further than the reciprocal distance of the pole from the $z=0$. The high-shelf is similar with the zero closer to the origin than the pole. The implementation posted by the OP was done with the zeros in the right half plane (in the s-plane) or equivalently outside the unit circle in the z-plane as I depict below (as the phase response posted by the OP goes to 180 degrees which tells us there are two zeros overall outside the unit circle). This results in a maximum phase response meaning maximum delay for the given magnitude response. It would likely be an improvement to instead reflect those zeros for each high-shelf and low-shelf section to be inside the unit circle as a minimum phase solution.
Other Approaches With Complex Conjugate Poles and Zeros
Independent control of bass or treble or both as the OP has done can be also be accomplished with a single biquad filter section with complex coefficients providing more control of the overall response. Given the entire phase is going to 180 degrees at Nyquist in all the OP's plots is the hint that suggests that these are all 2nd order systems of two poles and two zeros with all the poles inside the unit circle and all the zeros outside the unit circle.
The first plot looks like what you could achieve by placing the conjugate poles and zeros on the same angle (resonant frequency), with the radius of the pole inside the unit circle and the zero outside, symmetrically at less distance from the unit circle than the reciprocal of the pole. Add a constant gain and the magnitude response shown on the upper left plot (for frequency from $0$ to $2\pi$ or DC to the sampling rate, in contrast to the OP's plots that go from $0$ to $\pi$ or DC to Nyquist) comes close to that as posted by the OP:
By adjusting the relative pole and zero conjugate locations (radius and angle) with independent angles for the poles and zeros as well as independent radii, we can achieve different balances for the bass and treble with the single second order biquadratic filter, as above a static gain can be added to shift the entire magnitude plot up or down:
Bass Boost
Treble Boost
I detailed the relationship between the poles and zeros and the coefficients / implementations for the 2nd order biquad in this recent post.
Like what you see? I teach online DSP courses where I bring intuition together with the math involved for a deeper and more creative understanding of signal processing concepts. You can find the latest course listings at https://dsprelated.com/courses and https://ieeeboston.org/courses/ Next courses start at the end of Feb 2024, with an early sign-up discount until Feb 15!