differences between freqz and signal.scipy.freqz

[bnicenboim]

Hi, not sure if the bug is on gsignal side (but to be honest signal gives the same output) or scipy, but you might want to investigate what's going on:

Here they mostly agree (except on one frequency that I don't need to use):

system <- list(b = c(0.000731569700125585, 0, -0.00292627880050234, 0, 
           0.00438941820075351, 0, -0.00292627880050234, 0, 0.000731569700125585
), a = c(1, -7.04203950456817, 21.7283807572297, -38.3801847495478, 
         42.4586614063362, -30.1289923970756, 13.3939745471932, -3.41070179210893, 
         0.380901732541468))

scipy <- reticulate::import("scipy")
# I want the frequency response for a particular frequency:
gsignal::freqz(system$b, system$a, n = 5)$h[2] 
#> [1] -0.000976+0.1062529i
# one can get the frequency response for a particular frequency (in the units of fs, by default: 2*pi ) in this way with scipy.signal (with a non integer:)
scipy$signal$freqz(system$b, system$a, worN=.2 * pi)[[2]]
#> [1] -0.000976+0.1062529i
# This should be equivalent to gsignal:
scipy$signal$freqz(system$b, system$a, worN= 5L)[[2]]
#> [1]  3.255208e-06+0.000000e+00i -9.759891e-04+1.062529e-01i
#> [3]  3.328587e-03+2.658062e-03i  3.101129e-04+1.143374e-04i
#> [5]  1.305432e-05+2.075128e-06i
gsignal::freqz(system$b, system$a, n = 5)$h
#> [1]  0.000000e+00+0.000000e+00i -9.759891e-04+1.062529e-01i
#> [3]  3.328587e-03+2.658062e-03i  3.101129e-04+1.143374e-04i
#> [5]  1.305432e-05+2.075128e-06i
# but it doesn't quite agree on the first value:

But here they quite disagree:

# Now if I want to use "wierder" number, they stop gsignal and scipy.signal stop agreeing
f <- 0.000390625
scipy$signal$freqz(system$b, system$a, worN=f * pi)[[2]]
#> [1] -9.171918-0.454378i
gsignal::freqz(system$b, system$a, n = as.integer(1/f))$h[2] 
#> [1] -1.1738-239.1172i

^{Created on 2022-06-10 by the reprex package (v2.0.1)}

I'm planning to start using Sos rather ba, which supposed to be more stable, I'll report if I also find differences there.

[gjmvanboxtel]

Hi,

Thanks for your report - I'll investigate. Meanwhile, using Sos instead of b/a is a good idea. My guess is that you will also find discrepancies there though.

Best, Geert

[gjmvanboxtel]

Hi Bruno,

It turns out that there were two issues here. First, the core code of freqz in gsignal (and in signal) was from 2005. Since then, freqz in Octave (and probably also Matlab) has been updated, with a major update in 2009 which I think is related to the issue that you raised. I now changed the core code of gsignal::freqz to reflect the current Octave 1.4.2 version, along with a few other changes which I thought useful. I committed the changes to Github but not yet to CRAN.

This corrected most of the issue, and gsignal::freqz now gives mostly the same result as scipy.signal.freqz.

The other issue was related to numerical precision when calculating the frequency response for frequencies very close to zero. The issue is probably related to subtle differences in the polyval functions used on the different platforms (including Matlab/Octave). Consider this:

in R:

system <- list(b = c(0.000731569700125585, 0, -0.00292627880050234, 0, 
                     0.00438941820075351, 0, -0.00292627880050234, 0, 0.000731569700125585
), a = c(1, -7.04203950456817, 21.7283807572297, -38.3801847495478, 
         42.4586614063362, -30.1289923970756, 13.3939745471932, -3.41070179210893, 
         0.380901732541468))
f <- 0.000390625
scipy <- reticulate::import("scipy")
scipy$signal$freqz(system$b, system$a, worN = (0:10 * f * pi))[[2]]
freqz(system$b, system$a, n = (0:10 * f * pi))$h

Python

[1]  0.000003255208+0.0000000i -9.171918279952-0.4543777i  0.289642967128+0.9349324i  0.654130591478+0.7407555i  0.808561983987+0.5799574i
[6]  0.884203249641+0.4630334i  0.925470128964+0.3758969i  0.950378215676+0.3091221i  0.966323779360+0.2558589i  0.977055349168+0.2120209i
[11]  0.984428474866+0.1749158i

gsignal

[1] 0.00000330478+0.0000000i 5.62241522433-0.1330604i 0.28597835842+0.9374853i 0.65361770126+0.7449039i 0.80891682883+0.5810155i
 [6] 0.88410652455+0.4628956i 0.92548846048+0.3759203i 0.95043019490+0.3091576i 0.96634743621+0.2558696i 0.97705176697+0.2120185i
[11] 0.98444392905+0.1749208i

in Matlab/Octave

b = [0.000731569700125585, 0, -0.00292627880050234, 0, 0.00438941820075351, 0, -0.00292627880050234, 0, 0.000731569700125585];
a = [1, -7.04203950456817, 21.7283807572297, -38.3801847495478, 42.4586614063362, -30.1289923970756, 13.3939745471932, -3.41070179210893, 0.380901732541468];
f = 0.000390625;
h = freqz(b, a, [0:10] * pi * f)

h =

  Columns 1 through 4:

  0.0000 +        0i   -52.7670 - 151.0123i     0.3006 +   0.9272i     0.6540 +   0.7415i

Columns 5 through 8:

  0.8091 +   0.5814i     0.8842 +   0.4630i     0.9255 +   0.3759i     0.9505 +   0.3092i

Columns 9 through 11:

  0.9663 +   0.2559i     0.9770 +   0.2120i     0.9844 +   0.1749i

You see that the values start to converge quickly once you get further away from frequency 0. Because this is probably is related to issues that are buried deep within the mathematical functions used by these scripts (too mathematical for me anyway), I do not think that there is much more that I can do.

Hope this helps, Geert