ar1st0crat / NWaves

.NET DSP library with a lot of audio processing functions
MIT License
453 stars 71 forks source link

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

Open joa77 opened 1 year ago

joa77 commented 1 year ago

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 commented 1 year ago

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.: https://github.com/ar1st0crat/NWaves/blob/0728ca68ffac85450814203f4876c7c9be2df7cf/NWaves/Filters/Fda/DesignIirFilter.cs#L133

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