LMS自适应滤波算法的C++实现

破晓 · 发表于 2015-8-20 14:49:36

头文件：

/*
* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
*
* This program is free software; you can redistribute it and/or modify it
* under the terms of the GNU General Public License as published by the
* Free Software Foundation, either version 2 or any later version.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* This program is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
* more details. A copy of the GNU General Public License is available at:
* http://www.fsf.org/licensing/licenses
*/
/*****************************************************************************
* lms.h
*
* Least Mean Square Adaptive Filter.
*
* Least mean squares (LMS) algorithms are a class of adaptive filter used
* to mimic a desired filter by finding the filter coefficients that relate
* to producing the least mean squares of the error signal (difference
* between the desired and the actual signal). It is a stochastic gradient
* descent method in that the filter is only adapted based on the error at
* the current time.
*
* This file implement three types of the LMS algorithm: conventional LMS,
* algorithm, LMS-Newton algorhm and normalized LMS algorithm.
*
* Zhang Ming, 2010-10, Xi'an Jiaotong University.
*****************************************************************************/
#ifndef LMS_H
#define LMS_H
#include <vector.h>
#include <matrix.h>
namespace splab
{
template<typename Type>
Type lms( const Type&, const Type&, Vector<Type>&, const Type& );
template<typename Type>
Type lmsNewton( const Type&, const Type&, Vector<Type>&,
const Type&, const Type&, const Type& );
template<typename Type>
Type lmsNormalize( const Type&, const Type&, Vector<Type>&,
const Type&, const Type& );
#include <lms-impl.h>
}
// namespace splab
#endif
// LMS_H

复制代码

实现文件：

/*
* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
*
* This program is free software; you can redistribute it and/or modify it
* under the terms of the GNU General Public License as published by the
* Free Software Foundation, either version 2 or any later version.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* This program is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
* more details. A copy of the GNU General Public License is available at:
* http://www.fsf.org/licensing/licenses
*/
/*****************************************************************************
* lms-impl.h
*
* Implementation for LMS Adaptive Filter.
*
* Zhang Ming, 2010-10, Xi'an Jiaotong University.
*****************************************************************************/
/**
* The conventional LMS algorithm, which is sensitive to the scaling of its
* input "xn". The filter order p = wn.size(), where "wn" is the Weight
* Vector, and "mu" is the iterative setp size, for stability "mu" should
* belong to (0, Rr[Rxx]).
*/
template <typename Type>
Type lms( const Type &xk, const Type &dk, Vector<Type> &wn, const Type &mu )
{
int filterLen = wn.size();
static Vector<Type> xn(filterLen);
// update input signal
for( int i=filterLen; i>1; --i )
xn(i) = xn(i-1);
xn(1) = xk;
// get the output
Type yk = dotProd( wn, xn );
// update the Weight Vector
wn += 2*mu*(dk-yk) * xn;
return yk;
}
/**
* The LMS-Newton is a variant of the LMS algorithm which incorporate
* estimates of the second-order statistics of the environment signals is
* introduced. The objective of the algorithm is to avoid the slowconvergence
* of the LMS algorithm when the input signal is highly correlated. The
* improvement is achieved at the expense of an increased computational
* complexity.
*/
template <typename Type>
Type lmsNewton( const Type &xk, const Type &dk, Vector<Type> &wn,
const Type &mu, const Type &alpha, const Type &delta )
{
assert( 0 < alpha );
assert( alpha <= Type(0.1) );
int filterLen = wn.size();
Type beta = 1-alpha;
Vector<Type> vP(filterLen);
Vector<Type> vQ(filterLen);
static Vector<Type> xn(filterLen);
// initialize the Correlation Matrix's inverse
static Matrix<Type> invR = eye( filterLen, Type(1.0/delta) );
// update input signal
for( int i=filterLen; i>1; --i )
xn(i) = xn(i-1);
xn(1) = xk;
Type yk = dotProd(wn,xn);
// update the Correlation Matrix's inverse
vQ = invR * xn;
vP = vQ / (beta/alpha+dotProd(vQ,xn));
invR = (invR - multTr(vQ,vP)) / beta;
// update the Weight Vector
wn += 2*mu * (dk-yk) * (invR*xn);
return yk;
}
/**
* The conventional LMS is very hard to choose a learning rate "mu" that
* guarantees stability of the algorithm. The Normalised LMS is a variant
* of the LMS that solves this problem by normalising with the power of
* the input. For stability, the parameter "rho" should beong to (0,2),
* and "gamma" is a small number to prevent <Xn,Xn> == 0.
*/
template <typename Type>
Type lmsNormalize( const Type &xk, const Type &dk, Vector<Type> &wn,
const Type &rho, const Type &gamma )
{
assert( 0 < rho );
assert( rho < 2 );
int filterLen = wn.size();
static Vector<Type> sn(filterLen);
// update input signal
for( int i=filterLen; i>1; --i )
sn(i) = sn(i-1);
sn(1) = xk;
// get the output
Type yk = dotProd( wn, sn );
// update the Weight Vector
wn += rho*(dk-yk)/(gamma+dotProd(sn,sn)) * sn;
return yk;
}

复制代码

测试代码：

/*****************************************************************************
* lms_test.cpp
*
* LMS adaptive filter testing.
*
* Zhang Ming, 2010-10, Xi'an Jiaotong University.
*****************************************************************************/
#define BOUNDS_CHECK
#include <iostream>
#include <iomanip>
#include <lms.h>
using namespace std;
using namespace splab;
typedef double Type;
const int N = 50;
const int order = 1;
const int dispNumber = 10;
int main()
{
Vector<Type> dn(N), xn(N), yn(N), wn(order+1);
for( int k=0; k<N; ++k )
{
xn[k] = Type(cos(TWOPI*k/7));
dn[k] = Type(sin(TWOPI*k/7));
}
int start = max(0,N-dispNumber);
Type xnPow = dotProd(xn,xn)/N;
Type mu = Type( 0.1 / ((order+1)*xnPow) );
Type rho = Type(1.0), gamma = Type(1.0e-9), alpha = Type(0.08);
cout << "The last " << dispNumber << " iterations of Conventional-LMS:" << endl;
cout << "observed" << "\t" << "desired" << "\t\t" << "output" << "\t\t"
<< "adaptive filter" << endl << endl;
wn = 0;
for( int k=0; k<start; ++k )
yn[k] = lms( xn[k], dn[k], wn, mu );
for( int k=start; k<N; ++k )
{
yn[k] = lms( xn[k], dn[k], wn, mu );
cout << setiosflags(ios::fixed) << setprecision(4)
<< xn[k] << "\t\t" << dn[k] << "\t\t" << yn[k] << "\t\t";
for( int i=0; i<=order; ++i )
cout << wn[i] << "\t";
cout << endl;
}
cout << endl << endl;
cout << "The last " << dispNumber << " iterations of LMS-Newton:" << endl;
cout << "observed" << "\t" << "desired" << "\t\t" << "output" << "\t\t"
<< "adaptive filter" << endl << endl;
wn = 0;
for( int k=0; k<start; ++k )
yn[k] = lmsNewton( xn[k], dn[k], wn, mu, alpha, xnPow );
for( int k=start; k<N; ++k )
{
yn[k] = lmsNewton( xn[k], dn[k], wn, mu, alpha, xnPow );
cout << setiosflags(ios::fixed) << setprecision(4)
<< xn[k] << "\t\t" << dn[k] << "\t\t" << yn[k] << "\t\t";
for( int i=0; i<=order; ++i )
cout << wn[i] << "\t";
cout << endl;
}
cout << endl << endl;
cout << "The last " << dispNumber << " iterations of Normalized-LMS:" << endl;
cout << "observed" << "\t" << "desired" << "\t\t" << "output" << "\t\t"
<< "adaptive filter" << endl << endl;
wn = 0;
for( int k=0; k<start; ++k )
yn[k] = lmsNormalize( xn[k], dn[k], wn, rho, gamma );
for( int k=start; k<N; ++k )
{
yn[k] = lmsNormalize( xn[k], dn[k], wn, rho, gamma );
cout << setiosflags(ios::fixed) << setprecision(4)
<< xn[k] << "\t\t" << dn[k] << "\t\t" << yn[k] << "\t\t";
for( int i=0; i<=order; ++i )
cout << wn[i] << "\t";
cout << endl;
}
cout << endl << endl;
cout << "The theoretical optimal filter is:\t\t" << "-0.7972\t1.2788"
<< endl << endl;
return 0;
}