如何在Python中应用自适应过滤器 [英] How to apply an adaptive filter in Python

问题描述

我想在Python中应用自适应过滤器，但是在网上找不到有关如何实现这种算法的任何文档或示例.我熟悉使用scipy.signal工具箱设计静态"滤波器的情况，但是我不知道怎么做是设计自适应滤波器.

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.

为了澄清:我有一个记录的信号S，其中包含噪声.在此记录中，有一个我想访问的真实"函数，称为T.我也有一个T的估计.我想设计一个过滤器，以使过滤后的S和T之间的错误最小化.请注意，在这种情况下，静态过滤器没有用，因为我正在尝试过滤非平稳信号.

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.

推荐答案

这是Python中带有Numpy的基本LMS自适应过滤器.
欢迎发表评论，最欢迎使用测试用例.

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

""" 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()

这篇关于如何在Python中应用自适应过滤器的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！