How to apply an adaptive filter in Python
Question: I would like to apply an adaptive filter in Python, but can't find any
documentation or examples online of how to implement such an algorithm. I'm
familiar with designing "static" filters using the `scipy.signal` toolbox, but
what I don't know how to do is design an adaptive filter.

To clarify: I have a recorded signal `S` which contains noise. Within this
recording there is a "true" function that I would like to access, call this
`T`. I also have an estimate of `T`. I want to design a filter such that the
error between the filtered `S` and `T` is minimised. Note that in this case a
static filter is not useful, as I am trying to filter a nonstationary signal.


Answer: Here's a basic LMS adaptive filter in Python with Numpy.  
Comments are welcome, testcases most welcome.

![enter image description here](http://i.stack.imgur.com/ePu9a.png)

    
    
    """ lms.py: a simple python class for Least mean squares adaptive filter """
    
    from __future__ import division
    import numpy as np
    
    __version__ = "2013-08-29 aug denis"
    
    #...............................................................................
    class LMS:
        """ lms = LMS( Wt, damp=.5 )  Least mean squares adaptive filter
        in:
            Wt: initial weights, e.g. np.zeros( 33 )
            damp: a damping factor for swings in Wt
    
        # for t in range(1000):
    
        yest = lms.est( X, y [verbose=] )
        in: X: a vector of the same length as Wt
            y: signal + noise, a scalar
            optional verbose > 0: prints a line like "LMS: yest y c"
        out: yest = Wt.dot( X )
            lms.Wt updated
    
        How it works:
        on each call of est( X, y ) / each timestep,
        increment Wt with a multiple of this X:
            Wt += c X
        What c would give error 0 for *this* X, y ?
    
            y = (Wt + c X) . X
            =>
            c = (y  -  Wt . X)
                --------------
                   X . X
    
        Swings in Wt are damped a bit with a damping factor a.k.a. mu in 0 .. 1:
            Wt += damp * c * X
    
        Notes:
            X s are often cut from a long sequence of scalars, but can be anything:
            samples at different time scales, seconds minutes hours,
            or for images, cones in 2d or 3d x time.
    
    """
    
    # See also:
    #     http://en.wikipedia.org/wiki/Least_mean_squares_filter
    #     Mahmood et al. Tuning-free step-size adaptation, 2012, 4p
    # todo: y vec, X (Wtlen,ylen)
    
    #...............................................................................
        def __init__( self, Wt, damp=.5 ):
            self.Wt = np.squeeze( getattr( Wt, "A", Wt ))  # matrix -> array
            self.damp = damp
    
        def est( self, X, y, verbose=0 ):
            X = np.squeeze( getattr( X, "A", X ))
            yest = self.Wt.dot(X)
            c = (y - yest) / X.dot(X)
                # clip to cmax ?
            self.Wt += self.damp * c * X
            if verbose:
                print "LMS: yest %-6.3g   y %-6.3g   err %-5.2g   c %.2g" % (
                    yest, y, yest - y, c )
            return yest
    
    #...............................................................................
    if __name__ == "__main__":
        import sys
    
        filterlen = 10
        damp = .1
        nx = 500
        f1 = 40  # chirp
        noise = .05 * 2  # * swing
        plot = 0
        seed = 0
    
        exec( "\n".join( sys.argv[1:] ))  # run this.py n= ...  from sh or ipython
        np.set_printoptions( 2, threshold=100, edgeitems=10, linewidth=80, suppress=True )
        np.random.seed(seed)
    
        def chirp( n, f0=2, f1=40, t1=1 ):  # <-- your test function here
            # from $scipy/signal/waveforms.py
            t = np.arange( n + 0. ) / n * t1
            return np.sin( 2*np.pi * f0 * (f1/f0)**t )
    
        Xlong = chirp( nx, f1=f1 )
        # Xlong = np.cos( 2*np.pi * freq * np.arange(nx) )
        if noise:
            Xlong += np.random.normal( scale=noise, size=nx )  # laplace ...
        Xlong *= 10
    
        print 80 * "-"
        title = "LMS  chirp  filterlen %d  nx %d  noise %.2g  damp %.2g " % (
            filterlen, nx, noise, damp )
        print title
        ys = []
        yests = []
    
    #...............................................................................
        lms = LMS( np.zeros(filterlen), damp=damp )
        for t in xrange( nx - filterlen ):
            X = Xlong[t:t+filterlen]
            y = Xlong[t+filterlen]  # predict
            yest = lms.est( X, y, verbose = (t % 10 == 0) )
            ys += [y]
            yests += [yest]
    
        y = np.array(ys)
        yest = np.array(yests)
        err = yest - y
        averr = "av %.2g += %.2g" % (err.mean(), err.std())
        print "LMS yest - y:", averr
        print "LMS weights:", lms.Wt
        if plot:
            import pylab as pl
            fig, ax = pl.subplots( nrows=2 )
            fig.set_size_inches( 12, 8 )
            fig.suptitle( title, fontsize=12 )
            ax[0].plot( y, color="orangered", label="y" )
            ax[0].plot( yest, label="yest" )
            ax[0].legend()
            ax[1].plot( err, label=averr )
            ax[1].legend()
            if plot >= 2:
                pl.savefig( "tmp.png" )
            pl.show()
    

