[SciPy-User] signal.firls

[Gregory Allen]

I was surprised to find that there was no signal.firls function in scipy for designing FIR filters using least-squares error minimization.

I wrote this based on a post in this list [http://mail.scipy.org/pipermail/scipy-user/2009-November/023101.html], added band weights, and fleshed it out in the style of signal.remez. It’s not as full-featured as the matlab version, but it solved a need for me. :)

It could be pasted at the bottom of signal/fir_filter_design.py. I include it here in case it is of any use to someone else, or to be included in scipy.

Thanks,
-Greg

---

import numpy as np
from scipy.special import sinc

def firls(numtaps, bands, desired, weight=None, Hz=1):
    """
    FIR filter design using least-squares error minimization.

    Calculate the filter coefficients for the finite impulse response
    (FIR) filter which has the best approximation to the desired
    frequency response described by `bands` and `desired` in the
    least squares sense.
    
    Parameters
    ----------
    numtaps : int
        The number of taps in the FIR filter.  `numtaps` must be odd.
    bands : array_like
        A monotonic sequence containing the band edges in Hz.
        All elements must be non-negative and less than half the sampling
        frequency as given by `Hz`.
    desired : array_like
        A sequence half the size of bands containing the desired gain
        in each of the specified bands.
    weight : array_like, optional
        A relative weighting to give to each band region. The length of
        `weight` has to be half the length of `bands`.
    Hz : scalar, optional
        The sampling frequency in Hz. Default is 1.

    Returns
    -------
    out : ndarray
        A rank-1 array containing the coefficients of the optimal
        (in a least squares sense) filter.

    Example
    -------
    We want to construct a filter with a passband at 0.2-0.3 Hz, and
    stop bands at 0-0.1 Hz and 0.4-0.5 Hz. Note that this means that the
    behavior in the frequency ranges between those bands is unspecified and
    may overshoot.
    >>> from scipy import signal
    >>> bpass = signal.firls(71, [0, 0.1, 0.2, 0.3, 0.4, 0.5], [0, 1, 0])
    >>> freq, response = signal.freqz(bpass)
    >>> ampl = np.abs(response)
    >>> import matplotlib.pyplot as plt
    >>> fig = plt.figure()
    >>> ax1 = fig.add_subplot(111)
    >>> ax1.semilogy(freq/(2*np.pi), ampl, 'b-')  # freq in Hz
    >>> plt.show()
    """
    if numtaps%2 == 0:
        raise ValueError("numtaps must be odd.")
    L = (numtaps-1)//2
    
    # normalize bands and make it 2 columns
    bands = np.asarray(bands).flatten()/Hz
    if len(bands)%2 == 1:
        raise ValueError("bands must contain frequency pairs.")
    bands = bands.reshape(-1,2)
    
    # check remaining params
    if len(bands) != len(desired):
        raise ValueError("desired must have one entry per band.")
    if weight is None:
        weight = np.ones_like(desired)

    # set up the linear matrix equation to be solved, Ax = b
    k = np.arange(L+1)[np.newaxis]
    m = k.T
    A,b = 0,0
    for i, (f0,f1) in enumerate(bands):
        Ai = f1 * (sinc(2*(m+k)*f1) + sinc(2*(m-k)*f1)) \
           - f0 * (sinc(2*(m+k)*f0) + sinc(2*(m-k)*f0))
        bi = desired[i] * (2*f1*sinc(2*m*f1) - 2*f0*sinc(2*m*f0))
        A += Ai * abs(weight[i]**2)
        b += bi * abs(weight[i]**2)
        
    # solve and return
    x = np.linalg.solve(A,b).squeeze()
    h = np.hstack((x[:0:-1]/2, x[0], x[1:]/2))
    return h

-------------- next part --------------
A non-text attachment was scrubbed...
Name: smime.p7s
Type: application/pkcs7-signature
Size: 4876 bytes
Desc: not available
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20150428/9166ca22/attachment.bin>