问题描述
我想在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中应用自适应过滤器的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！