scipy.signal工具箱设计“静态”滤波器，但我不知道如何设计自适应滤波器。S。在这个记录中，有一个“真实”的函数T。我还有一个T的估计值。我想设计一个过滤器，使得经过滤波后的S和T之间的误差最小化。请注意，在这种情况下，静态滤波器是无用的，因为我正在尝试过滤一个非平稳信号。
""" 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()
y，它只是一个标量...但是假设我想使用自适应滤波器从一个函数中提取信号，并使用信号的估计值。我不清楚如何在你的脚本中实现这一点。 - allhands你应该检查Padasip库 https://pypi.python.org/pypi/padasip
这个库中实现了多个自适应滤波器。你可以直接使用它，或者研究一下那里的滤波器是如何创建的（源代码在Github上 https://github.com/matousc89/padasip）
如果你想测试任何已实现的自适应滤波器在你的情况下的适用性，请看一下这个教程Padasip自适应滤波器基础知识 - 噪声消除、系统识别和信号预测