[SciPy-User] [SciPy-user] least-square filter design
Sun Nov 1 14:26:55 CST 2009
Neal Becker wrote:
> Tom K. wrote:
>> Neal Becker wrote:
>>> Anyone have code for least square (minimum mean square error) FIR filter
> I'm looking for something like this:
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
h[L+1:] = a[1:]/2.
def plot_response(h, name):
H = np.fft.fft(h, 2000)
f = np.arange(2000)/2000.
setp(gca(), xlim=(0, .5))
if __name__ == '__main__':
h = firls(31, [(0, .2), (.3, .5)])
from matplotlib.pyplot import *
h = firls(51, [(0, .25), (.35, .5)], [0, 1])
h = firls(51, [(0, .1), (.2, .3), (.4, .5)], [0, 1, 0])
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.
More information about the SciPy-User