TF zeros, poles and gain from Butterworth are different to scipy z, p, k

[ToGoOrNotToGo]

I am comparing NWaves results with matlab and scipy results. TF zeros, poles and gain from Butterworth are different to scipy z, p, k

var bw = new Filters.Butterworth.BandPassFilter(lower / fs, upper / fs, 6); var (z, p, k) = (bw.Tf.Zeros, bw.Tf.Poles, bw.Tf.Gain);

z, p, k = butter( N=6, Wn=np.array([lower, upper]) / (fs / 2), btype='bandpass', analog=False, output='zpk')

The results from z and p are reverse and so gain is different. Is it okay anyway?

[ar1st0crat]

Ordering of zeros and poles in output vectors doesn't matter (it's just a set of complex numbers). Gain should also be the same (or, at least, very close to scipy's result). If the gain is very small, there may be some insignificant difference, but it's caused by round-off errors. What input parameters did you specify (lower, upper, fs)?

[ToGoOrNotToGo]

Yes, scipy seems more precise. Maybe the calculation of the warped frequencies and the transformation into lowpass, bandpass, highpass, or bandstop are different?

Here some test code: `var fs = 48000d; var downsamplingFactor = 50; var fsd = fs / downsamplingFactor; var lower = 11.220184543019641;//calculated fractional-octave band var upper = 14.125375446227549;

var bw = new Filters.Butterworth.BandPassFilter(lower / fsd, upper / fsd, 6); var (z, p, k) = (bw.Tf.Zeros, bw.Tf.Poles, bw.Tf.Gain);

var tf = DesignFilter.TfToSos(bw.Tf); var sosFilter = new FilterChain(tf);`

[ar1st0crat]

The calculation is essentially the same (as far as I can judge from reading a sciPy codebase). But in your particular case the filter is very narrow, so round-off errors become an important factor. Actually, in sciPy the frequency response for your parameters is not good too:

check #33 and #30

and this

[ar1st0crat]

UPD: anyway, the filtering will be OK if you do it as below, using SOS:

var fs = 48000d;
var downsamplingFactor = 50;
var fsd = fs / downsamplingFactor;
var lower = 11.220184543019641;//calculated fractional-octave band
var upper = 14.125375446227549;

var filter = new Filters.Butterworth.BandPassFilter(lower / fsd, upper / fsd, 6);
var sos = DesignFilter.TfToSos(filter.Tf);
var sosFilter = new FilterChain(sos);

// 1) SOS will mitigate the numeric effects:
// 2) re-estimate gain (this isn't available in ver.0.9.4 yet, but it's in current codebase)

var gain = sosFilter.EstimateGain();
var filteredSignal = sosFilter.ApplyTo(_signal, gain);

[ToGoOrNotToGo]

Good to know, I'll try the estimate gain solution.

But at the moment I'm stuck with the scipy 'y = filtfilt(num, den, x)' equivalent ZiFilter.ZeroPhase. It returns big differences. I think because of the pad type. Default in scipy is odd padding. NWaves has only zero padding...

[ToGoOrNotToGo]

Oh ... I'm just beginning to understand. ZiFilter.ZeroPhase has an "implicit" odd padding, right? I'm trying to figure out where the big differences come from...

[ToGoOrNotToGo]

I found one issue. In TransferFunction.Zi the 'zi' array length should be size - 1: var zi = new double[size - 1];

[ToGoOrNotToGo]

An other issue seems to be the padLength in ZiFilter.ZeroPhase. It should be: padLength = 3 * Math.Max(_a.Length, _b.Length);

instead of: padLength = 3 * (Math.Max(_a.Length, _b.Length) - 1);

[ar1st0crat]

Yes, for some reason, sciPy filtfilt default parameter is not compliant with MATLAB/NWaves filtfilt.

In MATLAB/NWaves default padlen is 3 * (max(len(a), len(b)) - 1). If you need to achieve the same results in sciPy as in MATLAB/NWaves, specify this value explicitly. Example:

double[] num = {1, 0.2, 0.5, -0.6, 0.3};
double[] den = {1, -0.2, 0.4, 0.1};

var filter = new ZiFilter(num, den);

// ============== lfilter_zi(num, den): ==================

double[] zi = filter.Tf.Zi;

// ============== filtfilt(num, den, np.arange(20), padlen=3*(5-1)): ================== 

var input = new DiscreteSignal(1, Enumerable.Range(0, 20));

var output = filter.ZeroPhase(input);

In MATLAB:

y = filtfilt([1 0.2 0.5 -0.6 0.3], [1 -0.2 0.4 0.1], 0:19)

In sciPy:

num = [1, 0.2, 0.5, -0.6, 0.3]
den = [1, -0.2, 0.4, 0.1]

x = np.arange(20)

zi = lfilter_zi(num, den)
# [ 0.07692308, -0.33846154, -0.40769231,  0.3]
# in NWaves: { 0.07692308, -0.33846154, -0.40769231,  0.3, 0 }

y = filtfilt(num, den, x, padlen=12)        # 3* (max(len(num), len(den)) - 1)

PS. As for Zi length, it's not an issue. One trailing element is added for convenience (in internal computations). Simply ignore the last element.

[ToGoOrNotToGo]

One more comment on this: After many tests and some learning, I think that the inaccuracies arise because the gain (k) is not taken into account in the calculations with zeros and poles (z, p) in TransferFunction or DesignFilter.

Anyway ... NWaves is great for studying signal analysis, implementing math algorithms, OOP, and working with music audio. For general DSP and scientific use cases, further improvements are required (e.g. the z, p, k and the analog to digital filter transformations (lp to hp, bp, bs)).

Thank you for this library and your very good explanations! I hope you continue :-)

[ar1st0crat]

Thank you for your kind words!

The gain takes part in computations (actually it's always reflected in TF numerator). The reasons for inaccuracies in this particular case is the specifics of your particularly narrow filter (MATLAB, sciPy and NWaves - all of them give different results for this filter parameters due to numeric errors). Just check out other filters (with wider passband) - you'll see NWaves/MATLAB/sciPy will produce equivalent results.