optimization - How to improve filter quality at low fs

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)
//poles
double p1 = 3180e-6; // 1/(2*pi*50.05Hz)
double p2 = 75e-6;    // 2212Hz

//zeros
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.

Juha

Answer

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:

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.

EDIT:

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: