[SciPy-User] [SciPy-user] least-square filter design

[Tom K.]

Neal Becker wrote:
> 
> Tom K. wrote:
> 
>> 
>> 
>> Neal Becker wrote:
>>> 
>>> Anyone have code for least square (minimum mean square error) FIR filter
>>> design?
>>> 
> I'm looking for something like this: 
> http://www.mathworks.com/access/helpdesk/help/toolbox/filterdesign/ref/firls.html
> 
> 

Here's something that works for odd length, symmetric filters with constant
magnitude per band - it doesn't support everything that MathWorks' firls
supports (e.g. design of even length, differentiators, antisymmetric
filters, sloping bands) but hopefully this meets your need.

import numpy as np
from scipy.special import sinc

def firls(N, f, D=None):
    """Least-squares FIR filter.
    N -- filter length, must be odd
    f -- list of tuples of band edges
       Units of band edges are Hz with 0.5 Hz == Nyquist
       and assumed 1 Hz sampling frequency
    D -- list of desired responses, one per band
    """
    if D is None:
        D = [1, 0]
    assert len(D) == len(f), "must have one desired response per band"
    assert N%2 == 1, 'filter length must be odd'
    L = (N-1)//2

    k = np.arange(L+1)
    k.shape = (1, L+1)
    j = k.T

    R = 0
    r = 0
    for i, (f0, f1) in enumerate(f):
        R += np.pi*f1*sinc(2*(j-k)*f1) - np.pi*f0*sinc(2*(j-k)*f0) + \
             np.pi*f1*sinc(2*(j+k)*f1) - np.pi*f0*sinc(2*(j+k)*f0)

        r += D[i]*(2*np.pi*f1*sinc(2*j*f1) - 2*np.pi*f0*sinc(2*j*f0))

    a = np.dot(np.linalg.inv(R), r)
    a.shape = (-1,)
    h = np.zeros(N)
    h[:L] = a[:0:-1]/2.
    h[L] = a[0]
    h[L+1:] = a[1:]/2.
    return h

def plot_response(h, name):
    H = np.fft.fft(h, 2000)
    f = np.arange(2000)/2000.
    figure()
    semilogy(f, abs(H))
    grid()
    setp(gca(), xlim=(0, .5))
    xlabel('frequency (Hz)')
    ylabel('magnitude')
    title(name)

if __name__ == '__main__':
    h = firls(31, [(0, .2), (.3, .5)])
    from matplotlib.pyplot import *
    plot_response(h, 'lowpass')

    h = firls(51, [(0, .25), (.35, .5)], [0, 1])
    plot_response(h, 'highpass')

    h = firls(51, [(0, .1), (.2, .3), (.4, .5)], [0, 1, 0])
    plot_response(h, 'bandpass')
    show()

-- 
View this message in context: http://old.nabble.com/least-square-filter-design-tp26083443p26154273.html
Sent from the Scipy-User mailing list archive at Nabble.com.