1
$\begingroup$

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
$\endgroup$
  • $\begingroup$ For audio application, an example of processing nth order IIR filter : sourceforge.net/p/equalizerapo/code/HEAD/tree/trunk/filters/… $\endgroup$ – Juha P Jul 11 at 4:24
  • $\begingroup$ What are your filter orders? $\endgroup$ – Matt L. Jul 11 at 8:08
  • $\begingroup$ @MattL, I ran the code listing and it gave this data: The order of the low-pass is 3 lp_nums = 5.8060e-04 -5.8058e-04 -5.8058e-04 5.8060e-04 lp_denoms = 1.00000 -2.99773 2.99549 -0.99776 The order of the band-pass is 3 mp_nums = 0.006641 -0.026527 0.033131 0.000000 -0.033131 0.026527 -0.006641 mp_denoms = 1.00000 -5.96991 14.85409 -19.71720 14.72613 -5.86749 0.97438 The order of the high-pass is 5 hp_nums = 0.90139 -4.50002 8.99311 -8.99311 4.50002 -0.90139 hp_denoms = 1.00000 -4.79342 9.18639 -8.79679 4.20814 -0.80431 $\endgroup$ – Juha P Jul 11 at 9:46
  • $\begingroup$ @JuhaP: thx for running the code! $\endgroup$ – Matt L. Jul 11 at 10:14
0
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks, this is the direction I was trying to go in (I think)! I'm not sure if I should post this as a separate question, but do I have to do anything special to apply this to a band-pass? Even though the order of my band-pass is the same as my low-pass, it has twice the number of coefficients. Do I just treat it as if it's twice the order when I apply direct form II structure? $\endgroup$ – silentTee Jul 12 at 18:28
  • $\begingroup$ @silentTee: You just use the bandpass coefficients and implement the equation given in my answer. $\endgroup$ – Matt L. Jul 12 at 18:41
1
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ 5 dB ripple in the passband is much too big for using the audio at the output, but if I understood correctly, the OP just compares the outputs of each filter using thresholding. The output is not used as audio. $\endgroup$ – Matt L. Jul 11 at 13:35
  • $\begingroup$ Fair point,I think I misread the application.Still, clean implementation is not a bad thing and especially if it doesn't cost extra :-) $\endgroup$ – Hilmar Jul 11 at 17:11
  • $\begingroup$ @MattL. That may not be very practical. These are the step responses for three 4th order Cauer with 0.05, 0.5, and 5 dB ripple, DC-scaled, and matching attenuations (-40.05, -40.5, -45 dB). True, their transition-band is optimized, so the delays will differ a bit, but the ripple is quite scary. $\endgroup$ – a concerned citizen Jul 13 at 15:41
  • $\begingroup$ @aconcernedcitizen: I don't know if its practical for the given application. My comment was triggered by the fact that the above answer is based on the assumption that the filter outputs are used as audio signals, which turned out not to be the case. $\endgroup$ – Matt L. Jul 13 at 16:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.