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

头文件：

/*

* 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

#include

namespace splab

{

template

Type lms( const Type&, const Type&, Vector&, const Type& );

template

Type lmsNewton( const Type&, const Type&, Vector&,

const Type&, const Type&, const Type& );

template

Type lmsNormalize( const Type&, const Type&, Vector&,

const Type&, const Type& );

#include

}

// 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

Type lms( const Type &xk, const Type &dk, Vector &wn, const Type &mu )

{

int filterLen = wn.size();

static Vector 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

Type lmsNewton( const Type &xk, const Type &dk, Vector &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 vP(filterLen);

Vector vQ(filterLen);

static Vector xn(filterLen);

// initialize the Correlation Matrix's inverse

static Matrix 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 == 0.

*/

template

Type lmsNormalize( const Type &xk, const Type &dk, Vector &wn,

const Type &rho, const Type &gamma )

{

assert( 0 < rho );

assert( rho < 2 );

int filterLen = wn.size();

static Vector 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

#include

#include

using namespace std;

using namespace splab;

typedef double Type;

const int N = 50;

const int order = 1;

const int dispNumber = 10;

int main()

{

Vector dn(N), xn(N), yn(N), wn(order+1);

for( int k=0; 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

yn[k] = lms( xn[k], dn[k], wn, mu );

for( int k=start; 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

yn[k] = lmsNewton( xn[k], dn[k], wn, mu, alpha, xnPow );

for( int k=start; 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

yn[k] = lmsNormalize( xn[k], dn[k], wn, rho, gamma );

for( int k=start; 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;

}

运行结果：

The last 10 iterations of Conventional-LMS:

observed desired output adaptive filter

-0.2225 -0.9749 -0.8216 -0.5683 1.0810

0.6235 -0.7818 -0.5949 -0.5912 1.0892

1.0000 -0.0000 0.0879 -0.6084 1.0784

0.6235 0.7818 0.6991 -0.5983 1.0947

-0.2225 0.9749 0.8156 -0.6053 1.1141

-0.9010 0.4339 0.2974 -0.6294 1.1082

-0.9010 -0.4339 -0.4314 -0.6289 1.1086

-0.2225 -0.9749 -0.8589 -0.6239 1.1291

0.6235 -0.7818 -0.6402 -0.6412 1.1353

1.0000 -0.0000 0.0667 -0.6543 1.1271

The last 10 iterations of LMS-Newton:

observed desired output adaptive filter

-0.2225 -0.9749 -0.9739 -0.7958 1.2779

0.6235 -0.7818 -0.7805 -0.7964 1.2783

1.0000 -0.0000 0.0006 -0.7966 1.2783

0.6235 0.7818 0.7816 -0.7966 1.2784

-0.2225 0.9749 0.9743 -0.7968 1.2787

-0.9010 0.4339 0.4334 -0.7971 1.2787

-0.9010 -0.4339 -0.4340 -0.7971 1.2787

-0.2225 -0.9749 -0.9747 -0.7971 1.2788

0.6235 -0.7818 -0.7816 -0.7972 1.2789

1.0000 -0.0000 0.0001 -0.7973 1.2789

The last 10 iterations of Normalized-LMS:

observed desired output adaptive filter

-0.2225 -0.9749 -0.9749 -0.7975 1.2790

0.6235 -0.7818 -0.7818 -0.7975 1.2790

1.0000 -0.0000 -0.0000 -0.7975 1.2790

0.6235 0.7818 0.7818 -0.7975 1.2790

-0.2225 0.9749 0.9749 -0.7975 1.2790

-0.9010 0.4339 0.4339 -0.7975 1.2790

-0.9010 -0.4339 -0.4339 -0.7975 1.2790

-0.2225 -0.9749 -0.9749 -0.7975 1.2790

0.6235 -0.7818 -0.7818 -0.7975 1.2790

1.0000 -0.0000 -0.0000 -0.7975 1.2790

The theoretical optimal filter is: -0.7972 1.2788

Process returned 0 (0x0) execution time : 0.109 s

Press any key to continue.