<commit_before>"""Filter design.
"""

import types
import warnings

import numpy
from numpy import atleast_1d, poly, polyval, roots, real, asarray, allclose, \
    resize, pi, absolute, logspace, r_, sqrt, tan, log10, arctan, arcsinh, \
    cos, exp, cosh, arccosh, ceil, conjugate, zeros, sinh
from numpy import mintypecode
from scipy import special, optimize
from scipy.misc import comb

class BadCoefficients(UserWarning):
    pass

abs = absolute

def findfreqs(num, den, N):
    ep = atleast_1d(roots(den))+0j
    tz = atleast_1d(roots(num))+0j

    if len(ep) == 0:
        ep = atleast_1d(-1000)+0j

    ez = r_['-1',numpy.compress(ep.imag >=0, ep,axis=-1), numpy.compress((abs(tz) < 1e5) & (tz.imag >=0),tz,axis=-1)]

    integ = abs(ez) < 1e-10
    hfreq = numpy.around(numpy.log10(numpy.max(3*abs(ez.real + integ)+1.5*ez.imag))+0.5)
    lfreq = numpy.around(numpy.log10(0.1*numpy.min(abs(real(ez+integ))+2*ez.imag))-0.5)

    w = logspace(lfreq, hfreq, N)
    return w

def freqs(b,a,worN=None,plot=None):
    """Compute frequency response of analog filter.

    Description:

       Given the numerator (b) and denominator (a) of a filter compute its
       frequency response.

               b[0]*(jw)**(nb-1) + b[1]*(jw)**(nb-2) + ... + b[nb-1]
       H(w) = --------------------------------------------------------
               a[0]*(jw)**(na-1) + a[1]*(jw)**(na-2) + ... + a[na-1]

    Inputs:

       b, a --- the numerator and denominator of a linear filter.
       worN --- If None, then compute at 200 frequencies around the interesting
                parts of the response curve (determined by pole-zero locations).
                If a single integer, the compute at that many frequencies.
                Otherwise, compute the response at frequencies given in worN.
    Outputs: (w,h)

       w -- The frequencies at which h was computed.
       h -- The frequency response.
    """
    if worN is None:
        w = findfreqs(b,a,200)
    elif isinstance(worN, types.IntType):
        N = worN
        w = findfreqs(b,a,N)
    else:
        w = worN
    w = atleast_1d(w)
    s = 1j*w
    h = polyval(b, s) / polyval(a, s)
    if not plot is None:
        plot(w, h)
    return w, h

def freqz(b, a=1, worN=None, whole=0, plot=None):
    """Compute frequency response of a digital filter.

    Description:

       Given the numerator (b) and denominator (a) of a digital filter compute
       its frequency response.

                  jw               -jw            -jmw
           jw  B(e)    b[0] + b[1]e + .... + b[m]e
        H(e) = ---- = ------------------------------------
                  jw               -jw            -jnw
               A(e)    a[0] + a[2]e + .... + a[n]e

    Inputs:

       b, a --- the numerator and denominator of a linear filter.
       worN --- If None, then compute at 512 frequencies around the unit circle.
                If a single integer, the compute at that many frequencies.
                Otherwise, compute the response at frequencies given in worN
       whole -- Normally, frequencies are computed from 0 to pi (upper-half of
                unit-circle.  If whole is non-zero compute frequencies from 0
                to 2*pi.

    Outputs: (w,h)

       w -- The frequencies at which h was computed.
       h -- The frequency response.

    """
    b, a = map(atleast_1d, (b,a))
    if whole:
        lastpoint = 2*pi
    else:
        lastpoint = pi
    if worN is None:
        N = 512
        w = numpy.arange(0,lastpoint,lastpoint/N)
    elif isinstance(worN, types.IntType):
        N = worN
        w = numpy.arange(0,lastpoint,lastpoint/N)
    else:
        w = worN
    w = atleast_1d(w)
    zm1 = exp(-1j*w)
    h = polyval(b[::-1], zm1) / polyval(a[::-1], zm1)
    if not plot is None:
        plot(w, h)
    return w, h

def tf2zpk(b,a):
    """Return zero, pole, gain (z,p,k) representation from a numerator,
    denominator representation of a linear filter.
    """
    b,a = normalize(b,a)
    b = (b+0.0) / a[0]
    a = (a+0.0) / a[0]
    k = b[0]
    b /= b[0]
    z = roots(b)
    p = roots(a)
    return z, p, k

def zpk2tf(z,p,k):
    """Return polynomial transfer function representation from zeros
    and poles

    Inputs:

      z, p --- sequences representing the zeros and poles.
      k --- system gain.

    Outputs: (b,a)

      b, a --- numerator and denominator polynomials.
    """
    z = atleast_1d(z)
    k = atleast_1d(k)
    if len(z.shape) > 1:
        temp = poly(z[0])
        b = zeros((z.shape[0], z.shape[1]+1), temp.dtype.char)
        if len(k) == 1:
            k = [k[0]]*z.shape[0]
        for i in range(z.shape[0]):
            b[i] = k[i] * poly(z[i])
    else:
        b = k * poly(z)
    a = poly(p)
    return b, a

def normalize(b,a):
    """Normalize polynomial representation of a transfer function.
    """
    b,a = map(atleast_1d,(b,a))
    if len(a.shape) != 1:
        raise ValueError, "Denominator polynomial must be rank-1 array."
    if len(b.shape) > 2:
        raise ValueError, "Numerator polynomial must be rank-1 or rank-2 array."
    if len(b.shape) == 1:
        b = asarray([b],b.dtype.char)
    while a[0] == 0.0 and len(a) > 1:
        a = a[1:]
    if allclose(b[:,0], 0, rtol=1e-14):
        warnings.warn("Badly conditionned filter coefficients (numerator): the "
                      "results may be meaningless", BadCoefficients)
        while allclose(b[:,0], 0, rtol=1e-14) and (b.shape[-1] > 1):
            b = b[:,1:]<commit_msg>Mention the warning in tf2zpk and zpk2tf.<commit_after>"""Filter design.
"""

import types
import warnings

import numpy
from numpy import atleast_1d, poly, polyval, roots, real, asarray, allclose, \
    resize, pi, absolute, logspace, r_, sqrt, tan, log10, arctan, arcsinh, \
    cos, exp, cosh, arccosh, ceil, conjugate, zeros, sinh
from numpy import mintypecode
from scipy import special, optimize
from scipy.misc import comb

class BadCoefficients(UserWarning):
    pass

abs = absolute

def findfreqs(num, den, N):
    ep = atleast_1d(roots(den))+0j
    tz = atleast_1d(roots(num))+0j

    if len(ep) == 0:
        ep = atleast_1d(-1000)+0j

    ez = r_['-1',numpy.compress(ep.imag >=0, ep,axis=-1), numpy.compress((abs(tz) < 1e5) & (tz.imag >=0),tz,axis=-1)]

    integ = abs(ez) < 1e-10
    hfreq = numpy.around(numpy.log10(numpy.max(3*abs(ez.real + integ)+1.5*ez.imag))+0.5)
    lfreq = numpy.around(numpy.log10(0.1*numpy.min(abs(real(ez+integ))+2*ez.imag))-0.5)

    w = logspace(lfreq, hfreq, N)
    return w

def freqs(b,a,worN=None,plot=None):
    """Compute frequency response of analog filter.

    Description:

       Given the numerator (b) and denominator (a) of a filter compute its
       frequency response.

               b[0]*(jw)**(nb-1) + b[1]*(jw)**(nb-2) + ... + b[nb-1]
       H(w) = --------------------------------------------------------
               a[0]*(jw)**(na-1) + a[1]*(jw)**(na-2) + ... + a[na-1]

    Inputs:

       b, a --- the numerator and denominator of a linear filter.
       worN --- If None, then compute at 200 frequencies around the interesting
                parts of the response curve (determined by pole-zero locations).
                If a single integer, the compute at that many frequencies.
                Otherwise, compute the response at frequencies given in worN.
    Outputs: (w,h)

       w -- The frequencies at which h was computed.
       h -- The frequency response.
    """
    if worN is None:
        w = findfreqs(b,a,200)
    elif isinstance(worN, types.IntType):
        N = worN
        w = findfreqs(b,a,N)
    else:
        w = worN
    w = atleast_1d(w)
    s = 1j*w
    h = polyval(b, s) / polyval(a, s)
    if not plot is None:
        plot(w, h)
    return w, h

def freqz(b, a=1, worN=None, whole=0, plot=None):
    """Compute frequency response of a digital filter.

    Description:

       Given the numerator (b) and denominator (a) of a digital filter compute
       its frequency response.

                  jw               -jw            -jmw
           jw  B(e)    b[0] + b[1]e + .... + b[m]e
        H(e) = ---- = ------------------------------------
                  jw               -jw            -jnw
               A(e)    a[0] + a[2]e + .... + a[n]e

    Inputs:

       b, a --- the numerator and denominator of a linear filter.
       worN --- If None, then compute at 512 frequencies around the unit circle.
                If a single integer, the compute at that many frequencies.
                Otherwise, compute the response at frequencies given in worN
       whole -- Normally, frequencies are computed from 0 to pi (upper-half of
                unit-circle.  If whole is non-zero compute frequencies from 0
                to 2*pi.

    Outputs: (w,h)

       w -- The frequencies at which h was computed.
       h -- The frequency response.

    """
    b, a = map(atleast_1d, (b,a))
    if whole:
        lastpoint = 2*pi
    else:
        lastpoint = pi
    if worN is None:
        N = 512
        w = numpy.arange(0,lastpoint,lastpoint/N)
    elif isinstance(worN, types.IntType):
        N = worN
        w = numpy.arange(0,lastpoint,lastpoint/N)
    else:
        w = worN
    w = atleast_1d(w)
    zm1 = exp(-1j*w)
    h = polyval(b[::-1], zm1) / polyval(a[::-1], zm1)
    if not plot is None:
        plot(w, h)
    return w, h

def tf2zpk(b,a):
    """Return zero, pole, gain (z,p,k) representation from a numerator,
    denominator representation of a linear filter.

    If some values of b are too close to 0, they are removed. In that case, a
    BadCoefficients warning is emitted.
    """
    b,a = normalize(b,a)
    b = (b+0.0) / a[0]
    a = (a+0.0) / a[0]
    k = b[0]
    b /= b[0]
    z = roots(b)
    p = roots(a)
    return z, p, k

def zpk2tf(z,p,k):
    """Return polynomial transfer function representation from zeros
    and poles

    Inputs:

      z, p --- sequences representing the zeros and poles.
      k --- system gain.

    Outputs: (b,a)

      b, a --- numerator and denominator polynomials.

    If some values of b are too close to 0, they are removed. In that case, a
    BadCoefficients warning is emitted.
    """
    z = atleast_1d(z)
    k = atleast_1d(k)
    if len(z.shape) > 1:
        temp = poly(z[0])
        b = zeros((z.shape[0], z.shape[1]+1), temp.dtype.char)
        if len(k) == 1:
            k = [k[0]]*z.shape[0]
        for i in range(z.shape[0]):
            b[i] = k[i] * poly(z[i])
    else:
        b = k * poly(z)
    a = poly(p)
    return b, a

def normalize(b,a):
    """Normalize polynomial representation of a transfer function.

    If values of b are too close to 0, they are removed. In that case, a
    BadCoefficients warning is emitted.
    """
    b,a = map(atleast_1d,(b,a))
    if len(a.shape) != 1:
        raise ValueError, "Denominator polynomial must be rank-1 array."
    if len(b.shape) > 2:
        raise ValueError, "Numerator polynomial must be rank-1 or rank-2 array."
    if len(b.shape) == 1:
        b = asarray([b],b.dtype.char)
    while a[0] == 0.0 and len(a) > 1:
        a = a[1:]
    if allclose(b[:,0], 0, rtol=1e-14):
        warnings.warn("Badly conditionned filter coefficients (numerator): the "
                      "results may be meaningless", BadCoefficients)
        while allclose(b[:,0], 0, rtol=1e-14) and (b.shape[-1] > 1):
            b = b[:,1:]
{'fraction_non_alphanumeric': 0.1476374494993852, 'fraction_numerical': 0.022922887756894432, 'mean_word_length': 1.958430761236685, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '640793', 'n_tokens_mistral': 4118, 'n_tokens_neox': 3880, 'n_words': 1193}