ar1st0crat/NWaves

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

joa77 opened this issue · 1 comments

joa77 commented

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?

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.:

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