How to use Octave ellip() coefficients for 3-band C++ crossover

Question

I have been working on making a 3-band audio crossover in C++, with the goal of activating/deactivating a different device depending on whether the low, mid, or high range of the input is above a threshold. I got as far as using GNU Octave to generate and graph out the low, band, and high pass filter coefficients using the Elliptic IIR method. I'm pleased with the coefficient characteristics, but this is my first DSP software project so I don't know how to use them to actually implement the filter program.

I looked at a bunch of articles on the topic online, and I saw people say using the Biquad Direct type 2 form is the most efficient approach for this situation. When I looked at the math on Wikipedia I understood the concepts somewhat but became confused as to how to implement it, partially because the order of my filters is greater than 2, and partially because comparing the formulas to other sources it seems like there might even be a mistake in the Wiki source??

I feel like I am incredibly close to understanding what to do, but I need that final nudge in the right direction.

I'm not sure if it is helpful at this stage, but here is what I have made in Octave so far:

function iir_elliptic_version
  clf
  close all

  pkg load signal;
  sample_freq = 192000;
  
  
  [lp_nums, lp_denoms] = Gen_Low_Pass_Filter_Coeffs(159, sample_freq)
  [mp_nums, mp_denoms] = Gen_Mid_Pass_Filter_Coeffs(160, 1999, sample_freq)
  [hp_nums, hp_denoms] = Gen_High_Pass_Filter_Coeffs(2000, sample_freq)
  
  f = 0:sample_freq/2;
  W = f*(2*pi/sample_freq);
  H1 = freqz(lp_nums, lp_denoms, W);
  H2 = freqz(mp_nums, mp_denoms, W);
  H3 = freqz(hp_nums, hp_denoms, W);
  plot(f, 20*log10(abs(H1)), f, 20*log10(abs(H2)), f, 20*log10(abs(H3)))
  axis([0 4000 -60 20])
  title('Filter Frequency Response of Elliptic IIR filters')
  xlabel('Frequency (Hz)')
  ylabel('Magnitude (dB)')
  grid on
endfunction

function [nums, denoms] = Gen_Low_Pass_Filter_Coeffs( passband_end,
                                              Fs)
  %This function generates the coefficients for a low pass filter for the 
  %low range given the passband end frequency, and the sampling rate.
  stopband_start = passband_end + 25;
  pass_ripple_dB = 5;
  atten_dB = 20;
  norm_pass_freq = passband_end/(Fs/2);
  norm_stop_freq = stopband_start/(Fs/2);
  [order, cutoff_radians] = ellipord(norm_pass_freq, norm_stop_freq, pass_ripple_dB, atten_dB);
  printf("The order of the low-pass is %d\n", order);
  [nums, denoms] = ellip(order, pass_ripple_dB, atten_dB, cutoff_radians);
endfunction

function [nums, denoms] = Gen_Mid_Pass_Filter_Coeffs( passband_start,
                                                      passband_end,
                                                      Fs)
  %This function generates the coefficients for a band pass filter for the 
  %mid range given the passband start and end frequencies, and the sampling rate.
  pass_freqs = [passband_start, passband_end];
  stop_freqs = [passband_start - 160, passband_end + 170];
  pass_ripple_dB = 5;
  atten_dB = 20;
  norm_pass_freqs = pass_freqs/(Fs/2);
  norm_stop_freqs = stop_freqs/(Fs/2);
  [order, cutoff_radians] = ellipord(norm_pass_freqs, norm_stop_freqs, pass_ripple_dB, atten_dB);
  printf("The order of the band-pass is %d\n", order);
  [nums, denoms] = ellip(order, pass_ripple_dB, atten_dB, cutoff_radians);
endfunction

function [nums, denoms] = Gen_High_Pass_Filter_Coeffs( passband_start,
                                              Fs)
  %This function generates the coefficients for a high pass filter for the 
  %high range given the passband start frequency, and the sampling rate.
  stopband_end = passband_start - 25;
  pass_ripple_dB = 5;
  atten_dB = 20;
  norm_pass_freq = passband_start/(Fs/2);
  norm_stop_freq = stopband_end/(Fs/2);
  [order, cutoff_radians] = ellipord(norm_pass_freq, norm_stop_freq, pass_ripple_dB, atten_dB);
  printf("The order of the high-pass is %d\n", order);
  [nums, denoms] = ellip(order, pass_ripple_dB, atten_dB, cutoff_radians, 'high');
endfunction

Answer 1

EDIT: In reaction to Hilmar's answer I would like to clarify that my answer below assumes that the filter outputs are not used as audio signals, but that the outputs are used for thresholding, measuring which band contains most energy. This is what I understood from reading the first paragraph of the question. If I misunderstood and if the filter outputs are used as audio signals, then some good advice can be found in Hilmar's answer.

---

From your specs and from running the code (thx to @JuhaP) it turns out that the resulting filters have very low orders (between $3$ and $5$). In this case, and if you can use double precision in your implementation, I wouldn't worry about splitting the transfer function into second-order sections. Just implement the filters in a direct-form II structure:

$$y[n]=b_0x[n]+b_1x[n-1]+\ldots + b_Nx[n-N]-a_1y[n-1]-\ldots\\ -a_Ny[n-N]\tag{1}$$

where $N$ is the filter order.

If for whatever reason you want to use second-order structures, you can use the Matlab/Octave command tf2sos(), which converts your a and b coefficients to the coefficients of second-order sections.

Answer 2

Not sure what your application requirements are but there seem to be some issues with your implementation.

1. Your pass band ripple of 5 dB is way too big. Something like 0.1dB would be much more typical
2. Elliptic filters are poor choice for cross overs. The create significant phase distortions in the cross over region that you can easily see if you sum your three bands. The sum is nowhere near "flat". Better choice would be a Butterworth filter for odd order (using the proper sign) or Linkwitz-Riley for even orders.
3. It also looks like you are trying to build a 3-band equalizer, not a cross-over. If that's the case the whole approach is wrong .
4. Your filter is WAY to steep which creates an enormous amount of time domain ringing. Your low-pass filter rings for more than 100 milliseconds.
5. Try to stay away from $[num,den]$ representation of filters. Use $[z,p,k]$ or $sos$ instead.
6. Implement filters as cascaded biquads. It's easy to do and there are dozens of code examples on the internet
7. Do NOT use "Direct Form II". That's pretty much the worst for audio processing. Use Direct Form I or Transposed Form II instead.

A quick explanation for the last three bullets. In audio processing, the poles of the filters are often extremely close to the unit circle. It's easy to get the transfer function coefficients from the poles and zeros, but the other direction is numerically difficult since it requires finding the roots of a polynomial where the roots are very close together.

In evaluating different filter topology it is useful to look at the transfer function from the input (or output) to the state variables. For Direct Form II (and transposed Form I) this is given by the pole-only transfer function, which can get extremely large (in audio). Even your seemingly harmless low pass at 159 Hz has a pole-only gain of a whopping 94 dB. That means that your state variables are 10000 larger than your input/output.