守望者--AIR技术交流

 找回密码
 立即注册

QQ登录

只需一步,快速开始

扫一扫,访问微社区

搜索
热搜: ANE FlasCC 炼金术
查看: 852|回复: 0

[音频分析] 各种类型RLS自适应滤波算法的C++实现

[复制链接]
  • TA的每日心情
    擦汗
    2018-4-10 15:18
  • 签到天数: 447 天

    [LV.9]以坛为家II

    1742

    主题

    2094

    帖子

    13万

    积分

    超级版主

    Rank: 18Rank: 18Rank: 18Rank: 18Rank: 18

    威望
    562
    贡献
    29
    金币
    52620
    钢镚
    1422

    开源英雄守望者

    发表于 2015-8-20 14:57:09 | 显示全部楼层 |阅读模式
    头文件:
    1. /*
    2. * Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
    3. *
    4. * This program is free software; you can redistribute it and/or modify it
    5. * under the terms of the GNU General Public License as published by the
    6. * Free Software Foundation, either version 2 or any later version.
    7. *
    8. * Redistribution and use in source and binary forms, with or without
    9. * modification, are permitted provided that the following conditions are met:
    10. *
    11. * 1. Redistributions of source code must retain the above copyright notice,
    12. *    this list of conditions and the following disclaimer.
    13. *
    14. * 2. Redistributions in binary form must reproduce the above copyright
    15. *    notice, this list of conditions and the following disclaimer in the
    16. *    documentation and/or other materials provided with the distribution.
    17. *
    18. * This program is distributed in the hope that it will be useful, but WITHOUT
    19. * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    20. * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
    21. * more details. A copy of the GNU General Public License is available at:
    22. * http://www.fsf.org/licensing/licenses
    23. */


    24. /*****************************************************************************
    25. *                                   rls.h
    26. *
    27. * Recursive Least Square Filter.
    28. *
    29. * The RLS adaptive filter recursively finds the filter coefficients that
    30. * minimize a weighted linear least squares cost function relating to the
    31. * input signals. This in contrast to other algorithms such as the LMS that
    32. * aim to reduce the mean square error.
    33. *
    34. * The input signal of RLS adaptive filter is considered deterministic,
    35. * while for the LMS and similar algorithm it is considered stochastic.
    36. * Compared to most of its competitors, the RLS exhibits extremely fast
    37. * convergence. However, this benefit comes at the cost of high computational
    38. * complexity, and potentially poor tracking performance when the filter to
    39. * be estimated changes.
    40. *
    41. * This file includes five types usually used RLS algorithms, they are:
    42. * conventional RLS (rls),       stabilised fast transversal RLS (sftrls),
    43. * lattice RLS (lrls),           error feedblck lattice RLS (eflrls),
    44. * QR based RLS (qrrls).
    45. *
    46. * Zhang Ming, 2010-10, Xi'an Jiaotong University.
    47. *****************************************************************************/


    48. #ifndef RLS_H
    49. #define RLS_H


    50. #include <vector.h>
    51. #include <matrix.h>


    52. namespace splab
    53. {

    54.     template<typename Type>
    55.     Type rls( const Type&, const Type&, Vector<Type>&,
    56.               const Type&, const Type& );

    57.     template<typename Type>
    58.     Type sftrls( const Type&, const Type&, Vector<Type>&,
    59.                  const Type&, const Type&, const string& );

    60.     template<typename Type>
    61.     Type lrls( const Type&, const Type&, Vector<Type>&,
    62.                const Type&, const Type&, const string& );

    63.     template<typename Type>
    64.     Type eflrls( const Type&, const Type&, Vector<Type>&,
    65.                  const Type&, const Type&, const string& );

    66.     template<typename Type>
    67.     Type qrrls( const Type&, const Type&, Vector<Type>&,
    68.                 const Type&, const string& );


    69.     #include <rls-impl.h>

    70. }
    71. // namespace splab


    72. #endif
    73. // RLS_H
    复制代码


    实现文件:
    1. /*
    2. * Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
    3. *
    4. * This program is free software; you can redistribute it and/or modify it
    5. * under the terms of the GNU General Public License as published by the
    6. * Free Software Foundation, either version 2 or any later version.
    7. *
    8. * Redistribution and use in source and binary forms, with or without
    9. * modification, are permitted provided that the following conditions are met:
    10. *
    11. * 1. Redistributions of source code must retain the above copyright notice,
    12. *    this list of conditions and the following disclaimer.
    13. *
    14. * 2. Redistributions in binary form must reproduce the above copyright
    15. *    notice, this list of conditions and the following disclaimer in the
    16. *    documentation and/or other materials provided with the distribution.
    17. *
    18. * This program is distributed in the hope that it will be useful, but WITHOUT
    19. * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    20. * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
    21. * more details. A copy of the GNU General Public License is available at:
    22. * http://www.fsf.org/licensing/licenses
    23. */


    24. /*****************************************************************************
    25. *                              rls-impl.h
    26. *
    27. * Implementation for RLS Filter.
    28. *
    29. * Zhang Ming, 2010-10, Xi'an Jiaotong University.
    30. *****************************************************************************/


    31. /**
    32. * The conventional RLS algorighm. The parameter "lambda" is the Forgetting
    33. * Factor, the smaller "lambda" is, the smaller contribution of previous samples.
    34. * This makes the filter more sensitive to recent samples, which means more
    35. * fluctuations in the filter co-efficients. Suggesting range is: [0.8, 1.0].
    36. * The parametet "delta" is the value to initialeze the inverse of the Auto-
    37. * Relation Matrix of input signal, which can be chosen as an estimation of
    38. * the input signal power.
    39. */
    40. template <typename Type>
    41. Type rls( const Type &xk, const Type &dk, Vector<Type> &wn,
    42.           const Type &lambda, const Type &delta )
    43. {
    44.     assert( Type(0.8) <= lambda );
    45.     assert( lambda <= Type(1.0) );

    46.     int filterLen = wn.size();
    47.     Vector<Type> vP(filterLen);
    48.     Vector<Type> vQ(filterLen);

    49.     static Vector<Type> xn(filterLen);
    50.     static Matrix<Type> invR = eye( filterLen, Type(1.0/delta) );

    51.     // updata input signal
    52.     for( int i=filterLen; i>1; --i )
    53.         xn(i) = xn(i-1);
    54.     xn(1) = xk;

    55.     // priori error
    56.     Type ak = dk - dotProd(wn,xn);

    57.     vQ = invR * xn;
    58.     vP = vQ / (lambda+dotProd(vQ,xn));

    59.     // updata Correlation-Matrix's inverse
    60.     invR = (invR - multTr(vQ,vP)) / lambda;

    61.     // update Weight Vector
    62.     wn += ak * vP;
    63.     //    wn += ak * (invR*xn);

    64.     return dotProd(wn,xn);
    65. }


    66. /**
    67. * Stabilized Fast Tranversal RLS
    68. */
    69. template <typename Type>
    70. Type sftrls( const Type &xk, const Type &dk, Vector<Type> &wn,
    71.              const Type &lambda, const Type &epsilon,
    72.              const string &training )
    73. {
    74.     int filterLen = wn.size(),
    75.         L = wn.size()-1;

    76.     assert( Type(1.0-1.0/(2*L+2)) <= lambda );
    77.     assert( lambda <= Type(1.0) );

    78.     static Vector<Type> xn(filterLen), xnPrev(filterLen);

    79.     const Type  k1 = Type(1.5),
    80.                 k2 = Type(2.5),
    81.                 k3 = Type(1.0);

    82.     // initializing for begin
    83.     Type    e, ep,
    84.             ef, efp,
    85.             eb1, eb2, ebp1, ebp2, ebp31, ebp32, ebp33;

    86.     static Type gamma = 1,
    87.                 xiBmin = epsilon,
    88.                 xiFminInv = 1/epsilon;

    89.     Vector<Type> phiExt(L+2);

    90.     static Vector<Type> phi(filterLen),
    91.                         wf(filterLen), wb(filterLen);

    92.     // updata input signal
    93.     xnPrev = xn;
    94.     for( int i=1; i<=L; ++i )
    95.         xn[i] = xnPrev[i-1];
    96.     xn[0] = xk;

    97.     if( training == "on" )
    98.     {
    99.         // forward prediction error
    100.         efp = xk - dotProd(wf,xnPrev);
    101.         ef  = gamma * efp;

    102.         phiExt[0] = efp * xiFminInv/lambda;
    103.         for( int i=0; i<filterLen; ++i )
    104.             phiExt[i+1] = phi[i] - phiExt[0]*wf[i];

    105.         // gamma1
    106.         gamma = 1 / ( 1/gamma + phiExt[0]*efp );

    107.         // forward minimum weighted least-squares error
    108.         xiFminInv = xiFminInv/lambda - gamma*phiExt[0]*phiExt[0];

    109.         // forward prediction coefficient vector
    110.         wf += ef * phi;

    111.         // backward prediction errors
    112.         ebp1 = lambda * xiBmin * phiExt[filterLen];
    113.         ebp2 = xnPrev[L] - dotProd(wb,xn);
    114.         ebp31 = (1-k1)*ebp1 + k1*ebp2;
    115.         ebp32 = (1-k2)*ebp1 + k2*ebp2;
    116.         ebp33 = (1-k3)*ebp1 + k3*ebp2;

    117.         // gamma2
    118.         gamma = 1 / ( 1/gamma - phiExt[filterLen]*ebp33 );

    119.         // backward prediction errors
    120.         eb1 = gamma * ebp31;
    121.         eb2 = gamma * ebp32;

    122.         // backward minimum weighted least-squares error
    123.         xiBmin = lambda*xiBmin + eb2*ebp2;

    124.         for( int i=0; i<filterLen; ++i )
    125.             phi[i] = phiExt[i] + phiExt[filterLen]*wb[i];

    126.         // forward prediction coefficient vector
    127.         wb += eb1 * phi;

    128.         // gamma3
    129.         gamma = 1 / ( 1 + dotProd(phi,xn) );

    130.         // Joint-Process Estimation
    131.         ep = dk - dotProd(wn,xn);
    132.         e  = gamma * ep;
    133.         wn += e * phi;
    134.     }

    135.     return dotProd(wn,xn);
    136. }


    137. /**
    138. * Lattice RLS
    139. */
    140. template <typename Type>
    141. Type lrls( const Type &xk, const Type &dk, Vector<Type> &vn,
    142.            const Type &lambda, const Type &epsilon,
    143.            const string &training )
    144. {
    145.     assert( Type(0.8) <= lambda );
    146.     assert( lambda <= Type(1.0) );

    147.     int filterLen = vn.size(),
    148.         L = filterLen-1;

    149.     // initializing for begin
    150.     Vector<Type>    gamma(filterLen),
    151.                     eb(filterLen),
    152.                     kb(L), kf(L),
    153.                     xiBmin(filterLen), xiFmin(filterLen);

    154.     static Vector<Type> delta(L), deltaD(filterLen),
    155.                         gammaOld(filterLen,Type(1.0)),
    156.                         ebOld(filterLen),
    157.                         xiBminOld(filterLen,epsilon),
    158.                         xiFminOld(filterLen,epsilon);

    159.     // initializing for Zeor Order
    160.     gamma[0] = 1;
    161.     xiBmin[0]= xk*xk + lambda*xiFminOld[0];
    162.     xiFmin[0] = xiBmin[0];

    163.     Type e = dk;
    164.     Type ef = xk;
    165.     eb[0] = xk;

    166.     for( int j=0; j<L; ++j )
    167.     {
    168.         // auxiliary parameters
    169.         delta[j] = lambda*delta[j] + ebOld[j]*ef/gammaOld[j];
    170.         gamma[j+1] = gamma[j] - eb[j]*eb[j]/xiBmin[j];

    171.         // reflection coefficients
    172.         kb[j] = delta[j] / xiFmin[j];
    173.         kf[j] = delta[j] / xiBminOld[j];

    174.         // prediction errors
    175.         eb[j+1] = ebOld[j] - kb[j]*ef;
    176.         ef -= kf[j]*ebOld[j];

    177.         // minimum least-squares
    178.         xiBmin[j+1] = xiBminOld[j] - delta[j]*kb[j];
    179.         xiFmin[j+1] = xiFmin[j] - delta[j]*kf[j];

    180.         // feedforward filtering
    181.         if( training == "on" )
    182.         {
    183.             deltaD[j] = lambda*deltaD[j] + e*eb[j]/gamma[j];
    184.             vn[j] = deltaD[j] / xiBmin[j];
    185.         }
    186.         e -= vn[j]*eb[j];
    187.     }

    188.     // last order feedforward filtering
    189.     if( training == "on" )
    190.     {
    191.         deltaD[L] = lambda*deltaD[L] + e*eb[L]/gamma[L];
    192.         vn[L] = deltaD[L] / xiBmin[L];
    193.     }
    194.     e -= vn[L]*eb[L];

    195.     // updated parameters
    196.     gammaOld = gamma;
    197.     ebOld = eb;
    198.     xiFminOld = xiFmin;
    199.     xiBminOld = xiBmin;

    200.     return dk-e;
    201. }


    202. /**
    203. * Error Feedback Lattice RLS
    204. */
    205. template <typename Type>
    206. Type eflrls( const Type &xk, const Type &dk, Vector<Type> &vn,
    207.              const Type &lambda, const Type &epsilon,
    208.              const string &training )
    209. {
    210.     assert( Type(0.8) <= lambda );
    211.     assert( lambda <= Type(1.0) );

    212.     int filterLen = vn.size(),
    213.         L = filterLen-1;

    214.     // initializing for begin
    215.     Vector<Type>    gamma(filterLen),
    216.                     eb(filterLen),
    217.                     xiBmin(filterLen), xiFmin(filterLen);

    218.     static Vector<Type> delta(L), deltaD(filterLen),
    219.                         gammaOld(filterLen,Type(1.0)),
    220.                         ebOld(filterLen),
    221.                         kb(L), kf(L),
    222.                         xiBminOld2(filterLen,epsilon),
    223.                         xiBminOld(filterLen,epsilon),
    224.                         xiFminOld(filterLen,epsilon);

    225.     // initializing for Zeor Order
    226.     gamma[0] = 1;
    227.     xiBmin[0]= xk*xk + lambda*xiFminOld[0];
    228.     xiFmin[0] = xiBmin[0];

    229.     Type tmp = 0;
    230.     Type e = dk;
    231.     Type ef = xk;
    232.     eb[0] = xk;

    233.     for( int j=0; j<L; ++j )
    234.     {
    235.         // auxiliary parameters
    236.         delta[j] = lambda*delta[j] + ebOld[j]*ef/gammaOld[j];
    237.         gamma[j+1] = gamma[j] - eb[j]*eb[j]/xiBmin[j];

    238.         // reflection coefficients
    239.         tmp = ebOld[j]*ef / gammaOld[j]/lambda;

    240.         kb[j] = gamma[j+1]/gammaOld[j] * ( kb[j] + tmp/xiFminOld[j] );
    241.             kf[j] = gammaOld[j+1]/gammaOld[j] * ( kf[j] + tmp/xiBminOld2[j] ) ;

    242.         // prediction errors
    243.         eb[j+1] = ebOld[j] - kb[j]*ef;
    244.         ef -= kf[j]*ebOld[j];

    245.         // minimum least-squares
    246.         xiBmin[j+1] = xiBminOld[j] - delta[j]*delta[j]/xiFmin[j];
    247.         xiFmin[j+1] = xiFmin[j] - delta[j]*delta[j]/xiBminOld[j];

    248.         // feedforward filtering
    249.         if( training == "on" )
    250.            vn[j] = gamma[j+1]/gamma[j] *
    251.                    ( vn[j] + e*eb[j]/(lambda*gamma[j]*xiBminOld[j]) );

    252.         e -= vn[j]*eb[j];
    253.     }

    254.     // last order feedforward filtering
    255.     if( training == "on" )
    256.         vn[L] = (gamma[L]-eb[L]*eb[L]/xiBmin[L])/gamma[L] *
    257.                 (vn[L]+e*eb[L]/(lambda*gamma[L]*xiBminOld[L]));
    258.     e -= vn[L]*eb[L];

    259.     // updated parameters
    260.     gammaOld = gamma;
    261.     ebOld = eb;
    262.     xiBminOld2 = xiBminOld;
    263.     xiBminOld = xiBmin;
    264.     xiFminOld = xiFmin;

    265.     return dk-e;
    266. }


    267. /**
    268. * QR-RLS
    269. */
    270. template <typename Type>
    271. Type qrrls( const Type &xk, const Type &dk, Vector<Type> &wn,
    272.             const Type &lambdaSqrt, const string &training )
    273. {
    274.     int filterLen = wn.size(),
    275.                 fL1 = filterLen+1;

    276.     // initializing for begin
    277.         static int k = 1;

    278.         Type        dp, ep,
    279.             gammap,
    280.                         c, cosTheta, sinTheta,
    281.                         tmp;

    282.     static Type  sx1;

    283.     Vector<Type> xp(filterLen);

    284.     static Vector<Type> xn(filterLen), dn(filterLen), dq2p(filterLen);

    285.         static Matrix<Type> Up(filterLen,filterLen);

    286.         // updata input signal
    287.         for( int i=filterLen; i>1; --i )
    288.             xn(i) = xn(i-1);
    289.         xn(1) = xk;

    290.     if( training == "on" )
    291.     {
    292.         // initializing for 0 to L iterations
    293.         if( k <= filterLen )
    294.         {
    295.             // updata Up
    296.             for( int i=filterLen; i>1; --i )
    297.                 for( int j=1; j<=filterLen; ++j )
    298.                     Up(i,j) = lambdaSqrt * Up(i-1,j);
    299.             for( int j=1; j<=filterLen; ++j )
    300.                     Up(1,j) = lambdaSqrt * xn(j);

    301.             dn(k) = dk;
    302.             for( int i=filterLen; i>1; --i )
    303.                 dq2p(i) = lambdaSqrt * dq2p(i-1);
    304.             dq2p(1) = lambdaSqrt * dn(k);

    305.             if( k == 1 )
    306.             {
    307.                 sx1 = xk;
    308.                 if( abs(sx1) > Type(1.0-6) )
    309.                     sx1 = Type(1.0);
    310.             }

    311.             // new Weight Vector
    312.             wn(1) = dn(1) / sx1;
    313.             if( k > 1 )
    314.                 for( int i=2; i<=k; ++i )
    315.                 {
    316.                     tmp = 0;
    317.                     for( int j=2; j<=i; ++j )
    318.                         tmp += xn(k-j+1)*wn(i-j+1);

    319.                     wn(i) = (-tmp+dn(i)) / sx1;
    320.                 }

    321.             ep = dk - dotProd(wn,xn);

    322.             k++;
    323.         }
    324.         else
    325.         {
    326.             xp = xn;
    327.             gammap = 1;
    328.             dp = dk;

    329.             // Givens rotation
    330.             for( int i=1; i<=filterLen; ++i )
    331.             {
    332.                 c = sqrt( Up(i,fL1-i)*Up(i,fL1-i) + xp(fL1-i)*xp(fL1-i) );
    333.                 cosTheta = Up(i,fL1-i) / c;
    334.                 sinTheta = xp(fL1-i) / c;
    335.                 gammap *= cosTheta;

    336.                 for( int j=1; j<=filterLen; ++j )
    337.                 {
    338.                     tmp = xp(j);
    339.                     xp(j) = cosTheta*tmp - sinTheta*Up(i,j);
    340.                     Up(i,j) = sinTheta*tmp + cosTheta*Up(i,j);
    341.                 }

    342.                 tmp = dp;
    343.                 dp = cosTheta*tmp - sinTheta*dq2p(i);
    344.                 dq2p(i) = sinTheta*tmp + cosTheta*dq2p(i);
    345.             }

    346.             ep = dp / gammap;

    347.             // new Weight Vector
    348.             wn(1) = dq2p(filterLen) / Up(filterLen,1);
    349.             for( int i=2; i<=filterLen; ++i )
    350.             {
    351.                 tmp = 0;
    352.                 for( int j=2; j<=i; ++j )
    353.                     tmp += Up(fL1-i,i-j+1) * wn(i-j+1);
    354.                 wn(i) = (-tmp+dq2p(fL1-i)) / Up(fL1-i,i);
    355.             }

    356.             // updating internal variables
    357.             Up *= lambdaSqrt;
    358.             dq2p *= lambdaSqrt;
    359.         }

    360.         return dk-ep;
    361.     }
    362.     else
    363.         return dotProd(wn,xn);
    364. }
    复制代码


    测试代码:
    1. /*****************************************************************************
    2. *                                   rls_test.cpp
    3. *
    4. * RLS adaptive filter testing.
    5. *
    6. * Zhang Ming, 2010-10, Xi'an Jiaotong University.
    7. *****************************************************************************/


    8. #define BOUNDS_CHECK

    9. #include <iostream>
    10. #include <iomanip>
    11. #include <convolution.h>
    12. #include <vectormath.h>
    13. #include <random.h>
    14. #include <rls.h>


    15. using namespace std;
    16. using namespace splab;


    17. typedef float   Type;
    18. const   int     N = 1000;
    19. const   int     orderRls = 1;
    20. const   int     orderLrls = 16;
    21. const   int     orderTrls = 12;
    22. const   int     orderQrrls = 8;
    23. const   int     sysLen = 8;
    24. const   int     dispNumber = 10;


    25. int main()
    26. {
    27.     int start = max(0,N-dispNumber);
    28.     Vector<Type>    dn(N), xn(N), yn(N), sn(N), rn(N), en(N),
    29.                     hn(sysLen+1), gn(orderLrls+1), wn(orderRls+1);
    30.     Type lambda, delta, eps;


    31.     cout << "/**************  Conventional RLS  <--->  Waveform Tracking \
    32. *************/" << endl << endl;
    33.     for( int k=0; k<N; ++k )
    34.     {
    35.         dn[k] = Type(sin(TWOPI*k/7));
    36.         xn[k] = Type(cos(TWOPI*k/7));
    37.     }
    38.     lambda = Type(0.99);
    39.     delta = dotProd(xn,xn)/N;

    40.     cout << "The last " << dispNumber << " iterations result:" << endl << endl;
    41.     cout << "observed" << "\t" << "desired" << "\t\t" << "output" << "\t\t"
    42.          << "adaptive filter" << endl << endl;
    43.     for( int k=0; k<start; ++k )
    44.         yn[k] = rls( xn[k], dn[k], wn, lambda, delta );
    45.     for( int k=start; k<N; ++k )
    46.     {
    47.         yn[k] = rls( xn[k], dn[k], wn, lambda, delta );
    48.         cout << setiosflags(ios::fixed) << setprecision(4)
    49.              << xn[k] << "\t\t" << dn[k] << "\t\t" << yn[k] << "\t\t";
    50.         for( int i=0; i<=orderRls; ++i )
    51.             cout << wn[i] << "\t";
    52.         cout << endl;
    53.     }
    54.     cout << endl << "The theoretical optimal filter is:\t\t" << "-0.7972\t1.2788"
    55.          << endl << endl << endl;


    56.     cout << "/**************  Lattice RLS  <--->  Channel Equalization \
    57. ***************/" << endl << endl;
    58.     for( int k=0; k<=sysLen; ++k )
    59.         hn[k] = Type( 0.1 * pow(0.5,k) );
    60.     dn = randn( 37, Type(0.0), Type(1.0), N );
    61.     xn = wkeep( conv(dn,hn), N, "left" );
    62.     lambda = Type(0.99), eps = Type(0.1);

    63.     wn.resize(orderLrls);
    64.     wn = Type(0.0);
    65.     for( int k=0; k<N; ++k )
    66. //        yn[k] = lrls( xn[k], dn[k], wn, lambda, eps, "on" );
    67.         yn[k] = eflrls( xn[k], dn[k], wn, lambda, eps, "on" );
    68.     Vector<Type> Delta(orderLrls+1);
    69.     Delta[(orderLrls+1)/2] = Type(1.0);
    70.     for( int k=0; k<=orderLrls; ++k )
    71. //        gn[k] = lrls( Delta[k], Type(0.0), wn, lambda, eps, "off" );
    72.         gn[k] = eflrls( Delta[k], Type(0.0), wn, lambda, eps, "off" );
    73.     cout << setiosflags(ios::fixed) << setprecision(4);
    74.     cout << "The original system:   " << hn << endl;
    75.     cout << "The inverse system:   " << gn << endl;
    76.     cout << "The cascade system:   " << conv( gn, hn ) << endl << endl;
    77. //

    78.     cout << "/************  Transversal RLS  <--->  System Identification \
    79. ************/" << endl << endl;
    80.     Vector<Type> sys(8);
    81.     sys[0] = Type(0.1);   sys[1] = Type(0.3);   sys[2] = Type(0.0);   sys[3] = Type(-0.2);
    82.     sys[4] = Type(-0.4);  sys[5] = Type(-0.7);  sys[6] = Type(-0.4);  sys[7] = Type(-0.2);
    83.     xn = randn( 37, Type(0.0), Type(1.0), N );
    84.     dn = wkeep( conv(xn,sys), N, "left" );
    85.     lambda = Type(0.99), eps = Type(1.0);
    86.     wn.resize(orderTrls);
    87.     wn = Type(0.0);

    88.     for( int k=0; k<start; ++k )
    89.         yn[k] = sftrls( xn[k], dn[k], wn, lambda, eps, "on" );
    90.     cout << "The last " << dispNumber << " iterations result:" << endl << endl;
    91.     cout << "input signal" << "    " << "original system output" << "    "
    92.          << "identified system output" << endl << endl;
    93.     for( int k=start; k<N; ++k )
    94.     {
    95.         yn[k] = sftrls( xn[k], Type(0.0), wn, lambda, eps, "off" );
    96.         cout << setiosflags(ios::fixed) << setprecision(4)
    97.              << xn[k] << "\t\t\t" << dn[k] << "\t\t\t" << yn[k] << endl;
    98.     }
    99.     cout << endl << "The unit impulse response of original system:   " << sys << endl;
    100.     cout << "The unit impulse response of identified system:   " << wn << endl << endl;


    101.     cout << "/***************  QR Based RLS  <--->  Signal Enhancement \
    102. ***************/" << endl << endl;
    103.     sn = sin( linspace( Type(0.0), Type(4*TWOPI), N ) );
    104.     rn = randn( 37, Type(0.0), Type(1.0), N );
    105.     dn = sn + rn;
    106.     int delay = orderQrrls/2;
    107.     for( int i=0; i<N-delay; ++i )
    108.         xn[i] = rn[i+delay];
    109.     for( int i=N-delay; i<N; ++i )
    110.         xn[i] = 0;
    111.     Type lambdaSqrt = Type(1.0);
    112.     wn.resize(orderQrrls);
    113.     wn = Type(0.0);

    114.     for( int k=0; k<N; ++k )
    115.         yn[k] = qrrls( xn[k], dn[k], wn, lambdaSqrt, "on" );
    116.     en = dn - yn;
    117.     for( int i=0; i<N/2; ++i )
    118.     {
    119.         sn[i] = 0;
    120.         rn[i] = 0;
    121.         en[i] = 0;
    122.     }
    123.     cout << "The last " << dispNumber << " iterations result:" << endl << endl;
    124.     cout << "noised signal\t\t" << "enhanced signal\t\t" << "original signal"
    125.          << endl << endl;
    126.     for( int k=start; k<N; ++k )
    127.         cout << setiosflags(ios::fixed) << setprecision(4)
    128.              << dn[k] << "\t\t\t" << en[k] << "\t\t\t" << sn[k] << endl;
    129.     cout << endl << "The SNR before denoising is:   "
    130.          << 20*(log10(norm(sn))/log10(norm(rn))) << " dB" << endl;
    131.     cout << "The SNR after denoising is:    "
    132.          << 20*(log10(norm(en))/log10(norm(sn-en))) << " dB" << endl << endl;


    133.     return 0;
    134. }
    复制代码


    运行结果:
    1. /**************  Conventional RLS  <--->  Waveform Tracking *************/

    2. The last 10 iterations result:

    3. observed        desired         output          adaptive filter

    4. -0.9010         0.4339          0.4339          -0.7975 1.2790
    5. -0.9010         -0.4339         -0.4339         -0.7975 1.2790
    6. -0.2225         -0.9749         -0.9749         -0.7975 1.2790
    7. 0.6235          -0.7818         -0.7818         -0.7975 1.2790
    8. 1.0000          0.0000          0.0000          -0.7975 1.2790
    9. 0.6235          0.7818          0.7818          -0.7975 1.2790
    10. -0.2225         0.9749          0.9749          -0.7975 1.2790
    11. -0.9010         0.4339          0.4339          -0.7975 1.2790
    12. -0.9010         -0.4339         -0.4339         -0.7975 1.2790
    13. -0.2225         -0.9749         -0.9749         -0.7975 1.2790

    14. The theoretical optimal filter is:              -0.7972 1.2788


    15. /**************  Lattice RLS  <--->  Channel Equalization ***************/

    16. The original system:   size: 9 by 1
    17. 0.1000
    18. 0.0500
    19. 0.0250
    20. 0.0125
    21. 0.0063
    22. 0.0031
    23. 0.0016
    24. 0.0008
    25. 0.0004

    26. The inverse system:   size: 17 by 1
    27. 0.6038
    28. -0.0026
    29. -0.0018
    30. -0.0012
    31. 0.0008
    32. 0.0011
    33. 0.0011
    34. -0.0009
    35. 8.7646
    36. -5.0002
    37. 0.0000
    38. 0.0015
    39. -0.0007
    40. -0.0004
    41. 0.0014
    42. 0.0027
    43. 0.0048

    44. The cascade system:   size: 25 by 1
    45. 0.0604
    46. 0.0299
    47. 0.0148
    48. 0.0073
    49. 0.0037
    50. 0.0020
    51. 0.0011
    52. 0.0005
    53. 0.8767
    54. -0.0618
    55. -0.0309
    56. -0.0153
    57. -0.0077
    58. -0.0039
    59. -0.0018
    60. -0.0006
    61. 0.0002
    62. -0.0016
    63. 0.0002
    64. 0.0001
    65. 0.0000
    66. 0.0000
    67. 0.0000
    68. 0.0000
    69. 0.0000


    70. /************  Transversal RLS  <--->  System Identification ************/

    71. The last 10 iterations result:

    72. input signal    original system output    identified system output

    73. 1.5173                  -1.7206                 -1.7206
    74. -0.9741                 -1.3146                 -1.3146
    75. -1.4056                 -2.1204                 -2.1204
    76. -0.9319                 -2.3165                 -2.3165
    77. -0.5763                 -2.1748                 -2.1748
    78. 0.4665                  -1.2228                 -1.2228
    79. 0.5959                  0.7623                  0.7623
    80. 0.5354                  1.7904                  1.7904
    81. -0.4990                 1.6573                  1.6573
    82. -1.1519                 0.4866                  0.4866

    83. The unit impulse response of original system:   size: 8 by 1
    84. 0.1000
    85. 0.3000
    86. 0.0000
    87. -0.2000
    88. -0.4000
    89. -0.7000
    90. -0.4000
    91. -0.2000

    92. The unit impulse response of identified system:   size: 12 by 1
    93. 0.1000
    94. 0.3000
    95. 0.0000
    96. -0.2000
    97. -0.4000
    98. -0.7000
    99. -0.4000
    100. -0.2000
    101. 0.0000
    102. 0.0000
    103. -0.0000
    104. -0.0000


    105. /***************  QR Based RLS  <--->  Signal Enhancement  ***************/

    106. The last 10 iterations result:

    107. noised signal           enhanced signal         original signal

    108. 1.2928                  -0.2263                 -0.2245
    109. -1.1740                 -0.1998                 -0.1999
    110. -1.5808                 -0.1748                 -0.1752
    111. -1.0823                 -0.1469                 -0.1504
    112. -0.7017                 -0.1039                 -0.1255
    113. 0.3660                  -0.0745                 -0.1005
    114. 0.5205                  -0.0628                 -0.0754
    115. 0.4851                  -0.0439                 -0.0503
    116. -0.5242                 -0.0220                 -0.0252
    117. -1.1519                 0.0058                  0.0000

    118. The SNR before denoising is:   17.5084 dB
    119. The SNR after denoising is:    121.1680 dB


    120. Process returned 0 (0x0)   execution time : 0.187 s
    121. Press any key to continue.
    复制代码




    本文来自:http://my.oschina.net/zmjerry/blog/8531
    守望者AIR技术交流社区(www.airmyth.com)
    回复

    使用道具 举报

    您需要登录后才可以回帖 登录 | 立即注册

    本版积分规则

    
    关闭

    站长推荐上一条 /4 下一条

    QQ|手机版|Archiver|网站地图|小黑屋|守望者 ( 京ICP备14061876号

    GMT+8, 2024-3-29 08:46 , Processed in 0.053839 second(s), 33 queries .

    守望者AIR

    守望者AIR技术交流社区

    本站成立于 2014年12月31日

    快速回复 返回顶部 返回列表