LMS自适应滤波算法的C++实现
头文件：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 |
/* * 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 |
实现文件：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 |
/* * 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; } |
测试代码：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 |
/***************************************************************************** * 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; } |
运行结果：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 |
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. |
