BUG: signal.iirfilter of bandpass type is unstable when high corner exactly at Nyquist

[megies]

We have come across what seems to be a bug in scipy.signal.iirfilter (or some underlying routine), see obspy/obspy#1451.

When the high corner of a digital butterworth bandpass filter happens to be exactly at Nyquist frequency, i.e. iirfilter(filter_order, (x, 1.0), btype='band', analog=False, ftype='butter', ..), (which is not raising an exception as is done when using a high corner above Nyquist), the filter gives bad results (the exact looks of which depending on a combination of filter order and low corner).

To illustrate the issue here's a plot and the code to generate it. Please excuse the longish example code, but I think the plot gives the best first impression of the issue. The signal is a real world example of a seismogram recording (link to data at the bottom).

• the first three rows use a high corner frequency just a tiny tad below Nyquist and show the expected results
• the bottom three rows use exactly Nyquist as high corner frequency and show wrong results, depending on the combination of low corner frequency and filter order..
  • show a signal resembling the original waveform (see below in black) or..
  • show a signal resembling the expected filter result but with a much lower amplitude

import numpy as np
import matplotlib.pyplot as plt
import scipy
from scipy.signal import iirfilter, sosfilt

data_unfiltered = np.loadtxt("seismogram_100Hz.txt")
sampling_rate = 100
nyquist = sampling_rate / 2.0

order = (3, 4, 5, 6)
high_corners = (49.99999, 50)
low_corners = (6.0, 8.55, 8.59)

fig = plt.figure(figsize=(10, 10))
num_axes_x = len(order)
num_axes_y = len(high_corners) * len(low_corners) + 1
num_axes = num_axes_x * num_axes_y

def _get_filtered_data(low, high, order):
    low = low / nyquist
    high = high / nyquist
    sos = iirfilter(order, [low, high], btype='band', ftype='butter',
                    output='sos')
    return sosfilt(sos, data_unfiltered.copy())

i = 0
ax = None
for high in high_corners:
    for low in low_corners:
        for order_ in order:
            data = _get_filtered_data(low=low, high=high, order=order_)
            ax = fig.add_subplot(num_axes_y, num_axes_x, i+1, sharex=ax,
                                 sharey=ax)
            ax.plot(data)
            if i < len(order):
                ax.set_title("order " + str(order_))
            if i % len(order) == 0:
                ax.set_ylabel("{} - {} Hz".format(low, high), rotation=0,
                              ha="right")
            i += 1

ax = plt.subplot(num_axes_y, num_axes_x, i+1, sharex=ax, sharey=ax)
ax.plot(data_unfiltered, "k-")
ax.set_ylabel("unfiltered/original data", rotation=0, ha="right")

plt.draw()
plt.suptitle("scipy version {} -- all x/y scales equal, Nyquist={} Hz, "
             "signal length is {} samples ({} seconds)".format(
                scipy.__version__, sampling_rate / 2.0, len(data_unfiltered),
                len(data_unfiltered) / sampling_rate))
plt.tight_layout()
fig.subplots_adjust(hspace=0, wspace=0, top=0.9)
for ax in fig.axes:
    ax.set_yticks([])
    ax.set_xticks([])
plt.show()

Data used in the example: seismogram_100Hz.txt

I don't think it matters but I'm on Linux 64bit using Anaconda, Python 2.7:

basemap                   1.0.8.dev0          np110py27_1    conda-forge
blas                      1.1                    openblas    conda-forge
cairo                     1.12.18                       1    conda-forge
curl                      7.45.0                        0    defaults
cycler                    0.10.0                   py27_0    defaults
decorator                 4.0.9                    py27_0    defaults
flake8                    2.5.1                    py27_0    defaults
fontconfig                2.11.1                        2    conda-forge
freetype                  2.6.3                         0    conda-forge
funcsigs                  0.4                      py27_0    http://repo.continuum.io/pkgs/free/linux-64/funcsigs-0.4-py27_0.tar.bz2
functools32               3.2.3.2                  py27_1    conda-forge
future                    0.14.3                   py27_0    http://repo.continuum.io/pkgs/free/linux-64/future-0.14.3-py27_0.tar.bz2
gdal                      2.0.0                    py27_1    defaults
geos                      3.4.2                         0    defaults
hdf4                      4.2.11                        0    defaults
hdf5                      1.8.15.1                      2    defaults
icu                       56.1                          2    conda-forge
ipython                   4.1.1                    py27_0    defaults
ipython-genutils          0.1.0                     <pip>
ipython_genutils          0.1.0                    py27_0    defaults
jpeg                      8d                            0    <unknown>
kealib                    1.4.5                         0    defaults
libgcc                    5.2.0                         0    defaults
libgdal                   2.0.0                         1    defaults
libgfortran               1.0                           0    defaults
libiconv                  1.14                          2    conda-forge
libnetcdf                 4.3.3.1                       3    defaults
libpng                    1.6.17                        0    defaults
libxml2                   2.9.2                         0    defaults
libxslt                   1.1.28                        0    <unknown>
lxml                      3.5.0                    py27_0    defaults
matplotlib                1.5.1               np110py27_1    conda-forge
mccabe                    0.3.1                    py27_0    defaults
mkl                       11.3.1                        0    defaults
mkl-rt                    11.1                         p0    defaults
mkl-service               1.1.0                   py27_p0    defaults
mock                      1.3.0                    py27_0    defaults
nose                      1.3.7                    py27_0    defaults
numexpr                   2.4.4               np110py27_0    defaults
numpy                     1.10.4          py27_blas_openblas_201  [blas_openblas]  conda-forge
obspy (/home/megies/git/obspy-maintenance_1.0.x) 1.0.1.post0+54.gcb3138c74f.obspy.fix.ascii.formats           <pip>
openblas                  0.2.18                        1    conda-forge
openssl                   1.0.2f                        0    defaults
pandas                    0.17.1              np110py27_0    defaults
path.py                   8.1.2                    py27_1    defaults
pbr                       1.3.0                    py27_0    defaults
pep8                      1.7.0                    py27_0    defaults
pexpect                   3.3                      py27_0    defaults
pickleshare               0.5                      py27_0    defaults
pip                       8.1.1                    py27_0    defaults
pixman                    0.32.6                        0    defaults
proj.4                    4.9.1                         0    defaults
proj4                     4.9.1                         0    defaults
py                        1.4.31                   py27_0    defaults
pycairo                   1.10.0                   py27_0    defaults
pyflakes                  1.0.0                    py27_0    defaults
pyparsing                 2.0.3                    py27_0    http://repo.continuum.io/pkgs/free/linux-64/pyparsing-2.0.3-py27_0.tar.bz2
pyproj                    1.9.4                    py27_1    defaults
pyqt                      4.11.4                   py27_1    defaults
pyshp                     1.2.3                    py27_0    conda-forge
pytest                    2.8.5                    py27_0    defaults
python                    2.7.11                        0    defaults
python-dateutil           2.4.2                    py27_0    defaults
pytz                      2015.7                   py27_0    defaults
qt                        4.8.7                         1    defaults
readline                  6.2                           2    <unknown>
requests                  2.9.1                    py27_0    defaults
scikit-learn              0.17.1          np110py27_blas_openblas_200  [blas_openblas]  conda-forge
scipy                     0.17.1          np110py27_blas_openblas_200  [blas_openblas]  conda-forge
setuptools                20.1.1                   py27_0    defaults
simplegeneric             0.8.1                    py27_0    defaults
sip                       4.16.9                   py27_0    defaults
six                       1.10.0                   py27_0    defaults
sqlalchemy                1.0.12                   py27_0    defaults
sqlite                    3.9.2                         0    defaults
tk                        8.5.18                        0    defaults
traitlets                 4.1.0                    py27_0    defaults
wheel                     0.29.0                   py27_0    defaults
xerces-c                  3.1.2                         0    defaults
zlib                      1.2.8                         0    http://repo.continuum.io/pkgs/free/linux-64/zlib-1.2.8-0.tar.bz2

[larsoner]

I suspect you might be hitting some boundary condition, where we probably shouldn't allow an upper limit equal to Nyquist (or we need to treat it differently to get it right). But I need to do some reading to figure that out...

In the meantime, you could look at the pole-zero plots of those filters, and also look at the frequency response with freqz to see how it varies.

[megies]

Figures, figures, figures:

scipy_filter_issue.py.txt

[e-q]

This may be the culprit: filter_design.py

if numpy.any(Wn < 0) or numpy.any(Wn > 1):
    raise ValueError("Digital filter critical frequencies must be 0 <= Wn <= 1")
fs = 2.0
warped = 2 * fs * tan(pi * Wn / fs)

[WarrenWeckesser]

The code that @e-q points out certainly is problematic if an element of Wn is 1. Theoretically, the warping doesn't work in this case, because Wn = 1 gets mapped to infinity, and the formulas used by the subsequent code are not designed to handle that case. In practice, because np.pi is not exactly π, tan(np.pi/2) is 16331239353195370.0. That's a big number, but still a poor approximation to infinity. :smile:

A quick fix is to simply not accept 1 in Wn. I haven't looked more closely to determine if there is a correct way to handle this edge case. A user can work around the problem by explicitly checking for 1, and designing a high-pass filter instead of a band-pass filter in that case.

[larsoner]

+1 for raising an error if Wn == 1, and letting the user change to a high-pass in that case. Won'

[larsoner]

...t make it for 0.18 but we should fix it before 0.19

[ev-br]

we should fix it before 0.19

Ping :-)

[larsoner]

Okay moved to 1.0 :)