This my 1st post here.

I'm using this impulse invariant (MZT) based method to calculate coefficients for biquad digital filter (for RIAA and non-RIAA de-emph.):

double a0, a1, a2, b0, b1, b2;
double fs = 44100;
//timeconstants (case RIAA):
// frequency -> time conversion 1/(2*pi*fc) (= R*C)
double p1 = 3180e-6; // 1/(2*pi*50.05Hz)
double p2 = 75e-6;    // 2212Hz
double z1 = 318e-6;  // 500.5Hz
double z2 = 0.0;       // 3.18e-6 for Neumann pole (50kHz)

double pole1= exp(-1.0/(fs*p1)); 
double pole2 = exp(-1.0/(fs*p2)); 
double zero1 = exp(-1.0/(fs*z1));
double zero2 = exp(-1.0/(fs*z2));

a0 = 1.0;   // = 1.0
a1 = -pole1 - pole2; // = -0.931176
a2 = pole1 * pole2; // = 0
b0 = 1.0; // = 1.0
b1 = -zero1 - zero2; // = -1.731986
b2 = zero1 * zero2; // = 0.733838

Known issue for MZT is magnitude error close Nyqvist at low sample rates - above filter results over -3dB error at 20kHz. I'm trying to find a way to fix this error. When using 96kHz sampling the error at 20kHz is around 0.5dB so the error (range 20Hz-20kHz) gets fixed by rising the sampling frequency high enough.

Not long ago Martin Visanec published an article "Matched Second Order Digital Filters" and he even kindly looked this case and tweaked his method suitable for above filter but the results he posted was not what I'm looking after (though, (without additional pole/zero pair) the result was less than ±0.3dB error in this above case).

Are there fitting methods one could consider to use for to improve the filter at low fs (for filtering system where coefficients are calculated in real time)?

Any thoughts and help with equations would be appreciated.


  • $\begingroup$ Which part of the specs is changing, such that you need to compute the filter coefficients in real time? I'm asking because I have a method for approximating analog filters by digital ones, but it is an iterative algorithm, no formula. $\endgroup$
    – Matt L.
    Mar 11, 2016 at 17:06
  • $\begingroup$ I'm working on this:jiiteepee.blogspot.com/2016/03/phonoeq-software-project.html. New filter is calculated when one of the three frequencies is changed. Also new coefficients is calculated if Neumann pole is selected (it's within the base coefficient as the code example shows). $\endgroup$
    – Juha P
    Mar 12, 2016 at 12:59
  • $\begingroup$ Actually, even the "Filter Gain" adjust changes the gain coefficients .... $\endgroup$
    – Juha P
    Mar 12, 2016 at 13:23

2 Answers 2


The impulse invariance method as well as the bilinear transform will give you relatively large errors close the Nyquist and there is nothing you can do about it. One thing you could do is use a filter design method that approximates the given analog response. I've tried this for the RIAA de-emphasis fitler with an (unpublished) algorithm I came up with some time ago, and the figure below shows the result:

enter image description here

The digital filter is a second-order IIR filter, the sampling frequency is $f_s=44100\text{ Hz}$, and the error is much less than 1dB over the whole frequency range up to Nyquist. The price you pay for it is 20 iterations of an approximation method, which solves a $5\times 5$ system of linear equations in each iteration.


If you add one more pole and one more zero (i.e. your discrete-time filter has 3 poles and 3 zeros), then the error goes down to less than 0.1 dB over the whole frequency range up to Nyquist:

enter image description here

  • $\begingroup$ @JuhaP: I haven't checked Martin Visanec's method, but note that I deliberately approximated up to Nyquist. You can reduce the error in the band up to 20kHz by excluding from the approximation the remaining band between 20kHz and Nyquist. Try FDLS and see what you get; it should be possible to use it with a single biquad. $\endgroup$
    – Matt L.
    Mar 12, 2016 at 13:23
  • $\begingroup$ @JuhaP: BTW, what was the reason you didn't want to use Visanec's method? $\endgroup$
    – Matt L.
    Mar 12, 2016 at 13:24
  • $\begingroup$ @JuhaP: There's no figure anymore in that link. I would be interested to see the result though. I've added a plot to my answer with one more pole and one more zero. The resulting error is now very small (< 0.1 dB over the whole frequency range). $\endgroup$
    – Matt L.
    Mar 12, 2016 at 17:27
  • $\begingroup$ @MattL. try i65.tinypic.com/mt48er.jpg $\endgroup$ Apr 13, 2016 at 10:05

Solved this by using binary search to find suitable position for leaky integrator.

  • $\begingroup$ This does not provide an answer to the question. To critique or request clarification from an author, leave a comment below their post. - From Review $\endgroup$
    – A_A
    Sep 18, 2016 at 6:49
  • $\begingroup$ @A_A This is the OP answering their own question. I do agree that it's not a particularly clear answer. If the OP can give more detail, that'd be great! $\endgroup$
    – Peter K.
    Sep 18, 2016 at 19:14
  • $\begingroup$ Thank you for letting me know. I do not remember typing this message but I do remember it poping up in my review queue where I tried to flag it as a "not-an-answer" as best as I could. I suppose the message was posted automatically (?) $\endgroup$
    – A_A
    Sep 18, 2016 at 19:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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