How to obtain digital SOS filter from analog zeros and poles?

[joa77]

I have a rather stupid question, but can't get my head around it: I have given zeros and poles (also a k value for gain correction) from an analog filter. I now want to create a SOS filter to filter audio signals with a given samplingrate. From my university time, I know that I have to use the bilinear transformation, to convert the analog poles&zeros to digital z domain, but yout BilinearTransform method, doesn't accept a samplerate parameter?

[ar1st0crat]

Your question is completely reasonable. The thing is, there's no particular method for the bilinear transform in NWaves. The BilinearTransform method simply performs complex division 1+z/1-z; this operation is at core of BT but it's not a bilinear transform per se (so the name is somewhat misleading). This method is called during IIR filter design (such as Butterworth, Chebyshev, etc.) where the "rest" of bilinear transform is implemented, e.g.:

NWaves/NWaves/Filters/Fda/DesignIirFilter.cs

Line 133 in 0728ca6

public static TransferFunction IirLpTf(double frequency, Complex[] poles, Complex[] zeros = null)

Note also that pre-warping is used: var warpedFreq = Math.Tan(Math.PI * frequency).

More details here

PS. The sampling rate is implicitly here because frequency parameter is the normalized frequency, i.e. freq / sampling_rate. Without pre-warping, according to MATLAB link above, the procedure would look something like this (didn't test it):

// input: zeros, poles, k

var dpoles = new Complex[PoleCount];
var dzeros = new Complex[ZeroCount];

// gain = real(k*prod(fs-z)./prod(fs-p));
var zprod = Complex.One;
var pprod = Complex.One;

// fs = 2*fs
sampling_rate *= 2;

for (var k = 0; k < poles.Length; k++)
{
    var p = poles[k] / sampling_rate;
    dpoles[k] = (1 + p) / (1 - p);

    pprod *= (sampling_rate - poles[k]);
}

for (var k = 0; k < zeros.Length; k++)
{
    var z = zeros[k] / sampling_rate;
    dzeros[k] = (1 + z) / (1 - z);

    zprod *= (sampling_rate - zeros[k]);
}

var gain = (k * zprod / pprod).Real;

var tf = new TransferFunction(dzeros, dpoles, gain);
// maybe also:
// tf.NormalizeAt(0);

TransferFunction[] sos = DesignFilter.TfToSos(tf);
var filter = new FilterChain(sos);

// apply filter
// ...