differences between ellip from gsignal and from scipy when the output is sos

[bnicenboim]

Hi, So I'm testing all my functions against python, and I'm seeing a different output between scipy ellip and gsignal ellip:

library(reticulate) 
library(gsignal)

signal <- import("scipy")$signal

# ba and g_ba are exactly the same
ba <- signal$ellip(13, 0.009, 80, 0.05, output='ba')
g_ba <- ellip(13, 0.009, 80, 0.05, output='Arma')

#these two are quite different:
sos <- signal$ellip(13, 0.009, 80, 0.05, output='sos')
g_sos <- gsignal::ellip(13, 0.009, 80, 0.05, output='Sos')

Python output (in R):

>  round(sos,2)
     [,1]  [,2] [,3] [,4]  [,5] [,6]
[1,]    0  0.00    0    1 -0.93 0.00
[2,]    1 -1.79    1    1 -1.87 0.88
[3,]    1 -1.93    1    1 -1.90 0.92
[4,]    1 -1.96    1    1 -1.93 0.95
[5,]    1 -1.97    1    1 -1.95 0.98
[6,]    1 -1.97    1    1 -1.96 0.99
[7,]    1 -1.97    1    1 -1.97 1.00

gsignal output:

$sos
     [,1]      [,2] [,3] [,4]       [,5]      [,6]
[1,]    1 -1.788482    1    1 -1.8700907 0.8782369
[2,]    1 -1.931183    1    1 -1.9034256 0.9174554
[3,]    1 -1.957998    1    1 -1.9331860 0.9524327
[4,]    1 -1.966718    1    1 -1.9527107 0.9752949
[5,]    1 -1.970118    1    1 -1.9642356 0.9886076
[6,]    1 -1.971358    1    1 -1.9715768 0.9967269
[7,]    1  1.000000    0    1 -0.9270557 0.0000000

$g
[1] 6.377097e-05

attr(,"class")
[1] "Sos"

Is this due to differences in conventions? (or is it a bug?) Is there a way to always reconstruct the python output from the gsignal output?

[gjmvanboxtel]

Hi Bruno,

There are two issues that might account for the differences.

First, scipy.signal and gsignal (and Matlab/Octave, BTW) differ in the order in which the sos are returned from the call to the ellip() function. Section order is supposed to be important for fixed-point applications, but not or less so for floating-point applications (Smith, J.O.S., 2007, see here), Secondly, scipy.signal incorporates the system gain in the first section, while gsignal (and Matlab/Octave) do not.

If in gsignal you do:

g_zpk <- gsignal::ellip(13, 0.009, 80, 0.05, output='Zpg')
g_sos2 <- gsignal::zp2sos(g_zpk$z, g_zpk$p, g_zpk$g, 'up')
g_sos2$sos[1, 1:3] <- g_sos2$sos[1, 1:3] * g_sos2$g
g_sos2$g <- 1
round(g_sos2$sos, 2)

you get

      [,1]   [,2]  [,3]  [,4]    [,5]  [,6]
[1,]    0  0.00    0    1 -0.93 0.00
[2,]    1 -1.97    1    1 -1.97 1.00
[3,]    1 -1.97    1    1 -1.96 0.99
[4,]    1 -1.97    1    1 -1.95 0.98
[5,]    1 -1.96    1    1 -1.93 0.95
[6,]    1 -1.93    1    1 -1.90 0.92
[7,]    1 -1.79    1    1 -1.87 0.88

These are the same values as the Python output, but the sections are in a different order. Gsignal returns the sections in the same order as Octave, and Matlab has a slightly different order. Although I have not looked into it in detail, my guess is that these differences are due to differences in the exact algorithm for pairing the poles to the nearest zeros.

For your case, the frequency responses and the zero-pole plots are the same, regardless of the order of the sections

freqz(Sos(sos, 1))
mtext("Frequency response scipy.signal", 3)
freqz(g_sos)
mtext("Frequency response gsignal", 3)
zplane(Sos(sos, 1), main = "Pole-zero plot scipy.signal")
zplane(g_sos, main = "Pole-zero plot gsignal")

However, the question of whether the difference in ordering of the sections leads to a difference in the quality of filtered signal for your specific case requires further investigation.

Best, Geert

[bnicenboim]

That is super useful! Thanks a lot! I'll check if there are major differences in the filtered signal then