BUG: signal.cheby1: wrong gain and cut-off frequency

[rcasperson]

Describe your issue.

At a low pass |G(0)| shall be 0 dB and and |G(f_cutoff)| shall be -3 dB according to the definition of the characteristics of a low pass filter.
At scipy.signal.cheby1 |G(0)| is either 0 dB or -rp depending on the order of the filter and |G(f_cutoff)| is -rp.

The same issue occurs at cheby1 high pass and band pass filters.

Reproducing Code Example

import numpy as np
from scipy.signal import cheby1,freqs
f = np.arange(0,10E3,10)
b,a = cheby1(4,5,f[100],'lowpass',True)
_,G = freqs(b,a,f)
for i in [0,100]: print(f"|G({f[i]:4.0f} Hz)| = {20*np.log10(np.abs(G[i])):2.0f} dB")

Error message

No error message but wrong values calculated.

SciPy/NumPy/Python version and system information

1.11.4 1.26.4 sys.version_info(major=3, minor=11, micro=7, releaselevel='final', serial=0)
lapack_armpl_info:
  NOT AVAILABLE
lapack_mkl_info:
    libraries = ['mkl_rt']
    library_dirs = ['C:/Users/rcaspers/AppData/Local/anaconda3\\Library\\lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['C:/Users/rcaspers/AppData/Local/anaconda3\\Library\\include']
lapack_opt_info:
    libraries = ['mkl_rt']
    library_dirs = ['C:/Users/rcaspers/AppData/Local/anaconda3\\Library\\lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['C:/Users/rcaspers/AppData/Local/anaconda3\\Library\\include']
blas_armpl_info:
  NOT AVAILABLE
blas_mkl_info:
    libraries = ['mkl_rt']
    library_dirs = ['C:/Users/rcaspers/AppData/Local/anaconda3\\Library\\lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['C:/Users/rcaspers/AppData/Local/anaconda3\\Library\\include']
blas_opt_info:
    libraries = ['mkl_rt']
    library_dirs = ['C:/Users/rcaspers/AppData/Local/anaconda3\\Library\\lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['C:/Users/rcaspers/AppData/Local/anaconda3\\Library\\include']
Supported SIMD extensions in this NumPy install:
    baseline = SSE,SSE2,SSE3
    found = SSSE3,SSE41,POPCNT,SSE42,AVX,F16C,FMA3,AVX2
    not found = AVX512F,AVX512CD,AVX512_SKX,AVX512_CLX,AVX512_CNL,AVX512_ICL

[j-bowhay]

@DietBru maybe you would be interested in taking a look?

[DietBru]

Depending on the application, the way low-pass filters are characterized varies quite a bit and is often inconsistent across the literature. The parametrization of cheby1() adheres to the established conventions (as described in Wikipedia): The allowed ripples are confined to the interval [-rp, 0] dB and the edge frequency of the pass-band is defined as the crossing of the level -rp dB. As documented even order filters have a gain of -rp dB and odd orders have unity gain.

Wikipedia also provides a formula relating the 3 dB cut-off frequency f_c to Wp:

N = 5  # order of filter
rp = 5  # width of ripple band in dB
eps = np.sqrt(10**(rp/10) - 1)
s = np.cosh(1 / np.cosh(1/eps) / N)
b, a = cheby1(N, rp, f_c * s, 'lowpass', analog=True)

To obtain unity gain at 0 Hz for even order filters, just multiply the output b by 10**(rp/20) or alternatively a by 10**(-rp/20).

For practical purposes, I can recommend the nice graphical filter design tool pyFDA which allows interactive twiddling of filter parameters to a required specification.

My opinion on modifying cheby1():

• I am -1 on adding an option to scale to unity gain at 0 Hz. IMO that use case is too specific and can be easily done by adding an additional line.
• I am -0 on adding an option for specifying the 3 dB cut-off frequencies instead of rp. In my experience filter design is usually a bit of twiddling where the cut-off frequency is just one constraint which needs to be balanced among many others.

[rcasperson]

@DietBru: Thanks for your fast reply.

I think, one ore two additional options like cutoff="3dB" and gain="unity" would be nice.

Your explanation and sample code should be added to the API reference at https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.cheby1.html. Than additional options are not necessary.

[lucascolley]

Your explanation and sample code should be added to the API reference

@rcasperson would you be interested in making a docs PR for this, for @DietBru to review?

[rcasperson]

My docs PR is an addition to the notes: The equiripple passband has N maxima or minima (for example, a 5th-order filter has 3 maxima and 2 minima). Consequently, the DC gain is unity for odd-order filters, or -rp dB for even-order filters. To obtain unity gain at 0 Hz for even-order filters, just multiply the output b by 10**(rp/20) or alternatively a by 10**(-rp/20). The edge frequency fe of the pass-band is defined as the crossing of the level -rp dB. If the cutoff frequency fc at a level of -3 dB is given, the required edge frequency at -rp can be calculated by fe = fc * s where s = np.cosh(1 / np.cosh(1/eps) / N) and eps = np.sqrt(10**(rp/10) - 1). Is it the same correction for low pass and high pass or is it fe = fc * s for a low pass and fe = fc / s for a high pass? Von: Lucas Colley ***@***.***> Gesendet: Donnerstag, 25. April 2024 11:16 An: scipy/scipy ***@***.***> Cc: Casperson, Ralf ***@***.***>; Mention ***@***.***> Betreff: Re: [scipy/scipy] BUG: signal.cheby1: wrong gain and cut-off frequency (Issue #20571) SICHERHEITSHINWEIS: Diese E-Mail wurde von außerhalb an die BAM gesendet. Bitte klicken Sie nicht auf Links oder öffnen Anhänge, bevor Sie nicht den Absender verifiziert haben und sicher sind, dass der Inhalt vertrauenswürdig ist. This message was sent from outside of BAM. Please do not click links or open attachments unless you recognize the sender and know the content is trustworthy. Your explanation and sample code should be added to the API reference @rcasperson <https://github.com/rcasperson> would you be interested in making a docs PR for this, for @DietBru <https://github.com/DietBru> to review? — Reply to this email directly, view it on GitHub <#20571 (comment)> , or unsubscribe <https://github.com/notifications/unsubscribe-auth/AE4VGEM3SYM3FVUKIJ2S3RLY7DCTFAVCNFSM6AAAAABGW3RZZWVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDANZWG4ZTGMJVGY> . You are receiving this because you were mentioned. <https://github.com/notifications/beacon/AE4VGEKPBD3QFGIIK5IDTHDY7DCTFA5CNFSM6AAAAABGW3RZZWWGG33NNVSW45C7OR4XAZNMJFZXG5LFINXW23LFNZ2KUY3PNVWWK3TUL5UWJTT3ZBXOI.gif> Message ID: ***@***.*** ***@***.***> >