守望者--AIR技术交流

 找回密码
 立即注册

QQ登录

只需一步,快速开始

扫一扫,访问微社区

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

[音频分析] Wiener滤波算法的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 15:09:39 | 显示全部楼层 |阅读模式
    头文件:
    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. *                                  levinson.h
    26. *
    27. * If the coefficient matrix of a linear equations is Toeplitz, then it can
    28. * be solved in a high computational efficiency way through Levinson-Durbin
    29. * algorithm. The subroutiones in this file will be used for solving Wiener
    30. * -Hopf equeations in Wiener filtring and Yule- Walker equations in
    31. * parametric spectrum estimation, respectively.
    32. *
    33. * Zhang Ming, 2010-11, Xi'an Jiaotong University.
    34. *****************************************************************************/


    35. #ifndef LEVINSON_H
    36. #define LEVINSON_H


    37. #include <vector.h>


    38. namespace splab
    39. {

    40.     template<typename Type>
    41.     Vector<Type> levinson( const Vector<Type>&, const Vector<Type>& );

    42.     template<typename Type>
    43.     Vector<Type> levinson( const Vector<Type>&, Type& );


    44.     #include <levinson-impl.h>

    45. }
    46. // namespace splab


    47. #endif
    48. // LEVINSON_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. *                                   wiener.h
    26. *
    27. * Wiener Filter.
    28. *
    29. * The goal of the Wiener filter is to filter out noise that has corrupted
    30. * a signal. It is based on a statistical approach.
    31. *
    32. * Wiener filters are characterized by the following:
    33. * Assumption:  signal and (additive) noise are stationary linear stochastic
    34. *              processes with known spectral characteristics or known
    35. *              autocorrelation and cross-correlation.
    36. * Requirement: the filter must be physically realizable/causal (this
    37. *              requirement can be dropped, resulting in a non-causal solution).
    38. * Performance: criterion: minimum mean-square error (MMSE).
    39. *
    40. * And The correlation matrix is a symmetric Toeplitz matrix, so we can use the
    41. * efficient algorithm, Levinson-Durbin algorithm, to solve the Wiener-Hopf
    42. * Equations.
    43. *
    44. * Zhang Ming, 2010-10, Xi'an Jiaotong University.
    45. *****************************************************************************/


    46. #ifndef WIENER_H
    47. #define WIENER_H


    48. #include <vector.h>
    49. #include <correlation.h>
    50. #include <levinson.h>


    51. namespace splab
    52. {

    53.     template<typename Type>
    54.     Vector<Type> wienerFilter( const Vector<Type>&,
    55.                                const Vector<Type>&, int );

    56.     template<typename Type>
    57.     Vector<Type> wienerPredictor( const Vector<Type>&, int );


    58.     #include <wiener-impl.h>

    59. }
    60. // namespace splab


    61. #endif
    62. // WIENER_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. *                               levinson-impl.h
    26. *
    27. * Implementationfor Levinson-Durbin alogrithm.
    28. *
    29. * Zhang Ming, 2010-11, Xi'an Jiaotong University.
    30. *****************************************************************************/


    31. /**
    32. * Levinson algorithm for solving Toeplitz equations.
    33. * t    : t(0), t(1), ..., t(n-1) of Toeplitz coefficient matrix
    34. * b    : constant vector
    35. */
    36. template <typename Type>
    37. Vector<Type> levinson( const Vector<Type> &t, const Vector<Type> &b )
    38. {
    39.     assert( t.size() == b.size() );

    40.     int n = t.size();
    41.         Type alpha, beta, q, c, omega;
    42.         Vector<Type> y(n), yy(n), x(n);

    43.         alpha = t[0];
    44.         if( abs(alpha) < EPS )
    45.         {
    46.                 cerr << "The matrix is ill-conditioned!" << endl;
    47.                 return x;
    48.         }
    49.         y[0] = 1;
    50.         x[0] = b[0] / alpha;

    51.         for( int k=1; k<n; ++k )
    52.         {
    53.                 q = 0;
    54.                 beta = 0;
    55.                 for( int j=0; j<k; ++j )
    56.                 {
    57.                         q += x[j] * t[k-j];
    58.                         beta += y[j] * t[j+1];
    59.                 }
    60.                 c = -beta / alpha;

    61.                 yy[0] = c * y[k-1];
    62.                 y[k] = y[k-1];
    63.                 for( int i=1; i<k; ++i )
    64.                         yy[i] = y[i-1] + c*y[k-i-1];
    65.                 yy[k] = y[k-1];

    66.                 alpha += c*beta;
    67.                 if( abs(alpha) < EPS )
    68.                 {
    69.                         cerr << "The matrix is ill-conditioned!" << endl;
    70.                         return x;
    71.                 }

    72.                 omega = (b[k]-q) / alpha;
    73.                 for( int i=0; i<k; ++i )
    74.                 {
    75.                         x[i] += omega*yy[i];
    76.                         y[i] = yy[i];
    77.                 }
    78.                 x[k] = omega*y[k];
    79.         }

    80.         return x;
    81. }


    82. /**
    83. * Levinson-Durbin algorithm for solving Youle-Walker equations.
    84. * rn       : r(0), r(1), ..., r(p)
    85. * sigma2   : the variance of exciting white noise
    86. */
    87. template <typename Type>
    88. Vector<Type> levinson( const Vector<Type> &rn, Type &sigma2 )
    89. {
    90.     int p = rn.size()-1;
    91.     Type tmp;
    92.     Vector<Type> ak(p+1), akPrev(p+1);

    93.     ak[0] = Type(1.0);
    94.     sigma2 = rn[0];
    95.     ak[1] = -rn[1]/sigma2;
    96.     sigma2 *= 1 - ak[1]*ak[1];

    97.     for( int k=2; k<=p; ++k )
    98.     {
    99.         tmp = 0;
    100.         for( int i=0; i<k; ++i )
    101.             tmp += ak[i]*rn[k-i];
    102.         ak[k] = -tmp/sigma2;

    103.         for( int i=1; i<k; ++i )
    104.             akPrev[i] = ak[i] + ak[k]*ak[k-i];
    105.         for( int i=1; i<k; ++i )
    106.             ak[i] = akPrev[i];

    107.         sigma2 *= 1 - ak[k]*ak[k];
    108.     }

    109.         return ak;
    110. }
    复制代码
    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. *                              wiener-impl.h
    26. *
    27. * Implementation for Wiener Filter.
    28. *
    29. * Zhang Ming, 2010-10, Xi'an Jiaotong University.
    30. *****************************************************************************/


    31. /**
    32. * By given the observed signal "xn" and desired signal "dn", this routine
    33. * return the coefficients for Wiener filter with "p" order.
    34. */
    35. template <typename Type>
    36. Vector<Type> wienerFilter( const Vector<Type> &xn,
    37.                            const Vector<Type> &dn, int p )
    38. {
    39.     int N = xn.size();
    40.     assert( dn.size() == N );

    41.     // auto-correlation and corss-correlation
    42.     Vector<Type> Rxx = fastCorr( xn, "unbiased" );
    43.     Vector<Type> Rdx = fastCorr( dn, xn, "unbiased" );

    44.     Vector<Type> tn(p+1), bn(p+1);
    45.     for( int i=0; i<=p; ++i )
    46.     {
    47.         tn[i] = Rxx[N-1+i];
    48.         bn[i] = Rdx[N-1+i];
    49.     }

    50.     // solving Wiener-Hopf equations
    51.     return levinson( tn, bn );
    52. }


    53. /**
    54. * One step Wiener predictor by using a "p" order filter. The input
    55. * signal "xn" should be much longer the "p" in order to have a more
    56. * precision estimation of the Correlation Matrix of "xn".
    57. */
    58. template <typename Type>
    59. Vector<Type> wienerPredictor( const Vector<Type> &xn, int p )
    60. {
    61.     int N = xn.size();

    62.     // auto-correlation and corss-correlation
    63.     Vector<Type> Rxx = fastCorr( xn, "unbiased" );
    64.     Vector<Type> tn(p+1), bn(p+1), predictor(p);

    65.     for( int i=0; i<=p; ++i )
    66.         tn[i] = Rxx[N-1+i];
    67.     bn(1) = Type(1.0);

    68.     // solving Yule-Walker equations
    69.     Vector<Type> tmp = levinson( tn, bn );

    70.     for( int i=1; i<=p; ++i )
    71.         predictor(i) = -tmp(i+1) / tmp(1);

    72.     return predictor;
    73. }
    复制代码


    测试代码:
    1. /*****************************************************************************
    2. *                                wiener_test.h
    3. *
    4. * Wiener 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 <wiener.h>
    12. #include <random.h>
    13. #include <vectormath.h>


    14. using namespace std;
    15. using namespace splab;


    16. typedef double  Type;
    17. const   int     N = 1024;
    18. const   int     M = 12;
    19. const   int     fOrder = 8;
    20. const   int     pOrder = 3;


    21. int main()
    22. {
    23.     Vector<Type> tn(N), dn(N), vn(N), xn(N), yn(N);
    24.     tn = linspace( Type(0.0), Type(2*TWOPI), N );
    25.     dn = sin(tn);
    26.     vn = randn( 37, Type(0.0), Type(1.0), N );
    27.     xn = dn + vn;

    28.     Vector<Type> hn = wienerFilter( xn, dn, fOrder );
    29.     cout << "Wiener filter hn:   " << hn << endl;
    30.     yn = wkeep( conv(xn,hn), N, "left" );
    31.     cout << "original relative error:   " << norm(dn-xn)/norm(dn) << endl;
    32.     cout << "filtered relative error:   " << norm(dn-yn)/norm(dn) << endl << endl;

    33.     Vector<Type> sn(M);
    34.     for( int i=0; i<M; ++i )
    35.         sn[i] = sin( Type(i*TWOPI/10) );
    36.     Vector<Type> pn = wienerPredictor( sn, pOrder );
    37.     Vector<Type> snPred = wkeep( conv(sn,pn), M, "left" );
    38.     cout << "Wiener predictor pn:   " << pn << endl;
    39.     cout << "original\tpredicted\terror" << endl;
    40.     Type realValue = sin( Type(M*TWOPI/10) );
    41.     cout << setiosflags(ios::fixed) << setprecision(4)
    42.          << realValue << "\t\t" << snPred[M-1] << "\t\t"
    43.          << abs(realValue-snPred[M-1]) << endl << endl;

    44.     return 0;
    45. }
    复制代码


    运行结果:
    1. Wiener filter hn:   size: 9 by 1
    2. 0.0903079
    3. 0.0858067
    4. 0.0813221
    5. 0.0870775
    6. 0.0968708
    7. 0.0888659
    8. 0.0844214
    9. 0.0901495
    10. 0.0965093

    11. original relative error:   1.43689
    12. filtered relative error:   0.416294

    13. Wiener predictor pn:   size: 3 by 1
    14. 0.606355
    15. 0.699992
    16. -1.03413

    17. original        predicted       error
    18. 0.9511          0.9643          0.0132


    19. Process returned 0 (0x0)   execution time : 0.094 s
    20. Press any key to continue.
    复制代码



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

    使用道具 举报

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

    本版积分规则

    
    关闭

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

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

    GMT+8, 2024-3-29 10:37 , Processed in 0.070413 second(s), 31 queries .

    守望者AIR

    守望者AIR技术交流社区

    本站成立于 2014年12月31日

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