守望者--AIR技术交流

 找回密码
 立即注册

QQ登录

只需一步,快速开始

扫一扫,访问微社区

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

[音频分析] 基于AR模型谱估计算法(Yule-Walker方法与Burg方法)的C++实现

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

    [LV.9]以坛为家II

    1742

    主题

    2094

    帖子

    13万

    积分

    超级版主

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

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

    开源英雄守望者

    发表于 2015-8-20 15:47:07 | 显示全部楼层 |阅读模式
    头文件:
    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. *                               parametricpse.h
    26. *
    27. * Parametric Power Spectrum Estimation Mothods.
    28. *
    29. * Techniques for spectrum estimation can generally be divided into parametric
    30. * (such as classical spectrum estimation) and non-parametric methods. The
    31. * parametric approaches assume that the underlying stationary stochastic
    32. * process has a certain structure which can be described using a small number
    33. * of parameters (for example, auto-regressive model). In these approaches,
    34. * the task is to estimate the parameters of the model that describes the
    35. * stochastic process.
    36. *
    37. * The widely used model is AR model, so this file provides three subroutines
    38. * to estimate the parameter of AR model, they are Yule-Walker method, Burg's
    39. * recursive mothod and forward-and-backward linear prediction least square
    40. * method.
    41. *
    42. * Zhang Ming, 2010-11, Xi'an Jiaotong University.
    43. *****************************************************************************/


    44. #ifndef PARAMETRICPSE_H
    45. #define PARAMETRICPSE_H


    46. #include <vector.h>
    47. #include <fft.h>
    48. #include <correlation.h>
    49. #include <levinson.h>
    50. #include <linequs1.h>


    51. namespace splab
    52. {

    53.     template<typename Type> Vector<Type> yulewalkerPSE( const Vector<Type>&,
    54.                                                         int, Type& );
    55.     template<typename Type> Vector<Type> burgPSE( const Vector<Type>&,
    56.                                                   int, Type& );
    57.     template<typename Type> Vector<Type> fblplsPSE( const Vector<Type>&,
    58.                                                     int, Type& );

    59.     template<typename Type> Vector<Type> armaPSD( const Vector<Type>&,
    60.                                                   const Vector<Type>&,
    61.                                                   const Type&, int );


    62.     #include <parametricpse-impl.h>

    63. }
    64. // namespace splab


    65. #endif
    66. // PARAMETRICPSE_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. *                            parametricpse-impl.h
    26. *
    27. * Implementation for parametric power spectrum estimatoin methods.
    28. *
    29. * Zhang Ming, 2010-11, Xi'an Jiaotong University.
    30. *****************************************************************************/


    31. /**
    32. * The Yule-Walker method for AR power spectral estimation.
    33. * xn       : input signal
    34. * p        : the AR model order
    35. * sigma2   : the variance of exciting white noise
    36. * return   : coefficients of AR model --- a(0), a(1), ..., a(p)
    37. */
    38. template <typename Type>
    39. Vector<Type> yulewalkerPSE( const Vector<Type> &xn, int p, Type &sigma2 )
    40. {
    41.     int N = xn.size();

    42.     assert( p <= N );

    43.     Vector<Type> rn(p+1);
    44.     for( int i=0; i<=p; ++i )
    45.         for( int k=0; k<N-i; ++k )
    46.             rn[i] += xn[k+i]*xn[k];
    47.     rn /= Type(N);

    48.     return levinson( rn, sigma2 );
    49. }


    50. /**
    51. * The Burg method for AR power spectral estimation.
    52. * xn       : input signal
    53. * p        : the AR model order
    54. * sigma2   : the variance of exciting white noise
    55. * return   : coefficients of AR model --- a(0), a(1), ..., a(p)
    56. */
    57. template <typename Type>
    58. Vector<Type> burgPSE( const Vector<Type> &xn, int p, Type &sigma2 )
    59. {
    60.     int N = xn.size();
    61.     Type numerator, denominator;
    62.     Vector<Type> ak(p+1), akPrev(p+1), ef(N), eb(N);

    63.     ak[0] = Type(1.0);
    64.     sigma2 = sum(xn*xn) / Type(N);

    65.     for( int i=1; i<N; ++i )
    66.     {
    67.         ef[i] = xn[i];
    68.         eb[i-1] = xn[i-1];
    69.     }

    70.     for( int k=1; k<=p; ++k )
    71.     {
    72.         numerator = 0;
    73.         denominator = 0;
    74.         for( int i=k; i<N; ++i )
    75.         {
    76.             numerator += ef[i]*eb[i-1];
    77.             denominator += ef[i]*ef[i] + eb[i-1]*eb[i-1];
    78.         }
    79.         ak[k] = -2*numerator/denominator;

    80.         for( int i=1; i<k; ++i )
    81.             akPrev[i] = ak[i] + ak[k]*ak[k-i];
    82.         for( int i=1; i<k; ++i )
    83.             ak[i] = akPrev[i];

    84.         sigma2 *= 1 - ak[k]*ak[k];

    85.         for( int i=N-1; i>k; --i )
    86.         {
    87.             ef[i] = ef[i] + ak[k]*eb[i-1];
    88.             eb[i-1] = eb[i-2] + ak[k]*ef[i-1];
    89.         }
    90.     }

    91.     return ak;
    92. }


    93. /**
    94. * Forward and backward linear prediction least square method for
    95. * AR power spectral estimation.
    96. * xn       : input signal
    97. * p        : the AR model order
    98. * sigma2   : the variance of exciting white noise
    99. * return   : coefficients of AR model --- a(0), a(1), ..., a(p)
    100. */
    101. template <typename Type>
    102. Vector<Type> fblplsPSE( const Vector<Type> &xn, int p, Type &sigma2 )
    103. {
    104.     int N = xn.size(),
    105.         M = 2*(N-p);

    106.     Vector<Type> u(p+1);
    107.     u[0] = Type(1.0);

    108.     Matrix<Type> X(M,p+1);
    109. //    for( int i=0; i<=p; ++i )
    110. //    {
    111. //        for( int j=0; j<N-p; ++j )
    112. //            X[j][i] = xn[i+j];
    113. //        for( int j=p; j<N; ++j )
    114. //            X[j][i] = xn[j-i];
    115. //    }
    116.     for( int i=0; i<N-p; ++i )
    117.         for( int j=0; j<=p; ++j )
    118.             X[i][j] = xn[i+j];
    119.     for( int i=p; i<N; ++i )
    120.         for( int j=0; j<=p; ++j )
    121.             X[i][j] = xn[i-j];

    122.     Matrix<Type> Rp = trMult(X,X) / Type(M);
    123.     Vector<Type> ak = luSolver( Rp, u );
    124.     sigma2 = 1/ak[0];
    125.     ak *= sigma2;

    126.     return ak;
    127. }


    128. /**
    129. * The Bartlett method of power spectral estimation.
    130. * ak       : AR coefficients
    131. * bk       : MA coefficients
    132. * sigma2   : the variance of exciting white noise
    133. * L        : the points number of PSD
    134. * return   : spectral density at L frequencies:
    135. *            w = 0, 2*pi/L, ..., 2*pi(L-1)/L
    136. */
    137. template <typename Type>
    138. Vector<Type> armaPSD( const Vector<Type> &ak, const Vector<Type> &bk,
    139.                       const Type &sigma2, int L )
    140. {
    141.     int p = ak.size()-1,
    142.         q = bk.size()-1;
    143.     Vector<Type> Xk(L);

    144.     Type zRe, zIm, aRe, aIm, bRe, bIm,
    145.          Xre, Xim,
    146.          re, im;
    147.         Type omega,
    148.          den, numRe, numIm;

    149.         for( int k=0; k<L; ++k )
    150.         {
    151.                 omega = Type(TWOPI*k/L);
    152.                 zRe = cos(-omega);
    153.                 zIm = sin(-omega);

    154.         // numerator
    155.                 bRe = 0;
    156.                 bIm = 0;
    157.                 for( int i=q; i>0; --i )
    158.                 {
    159.                         re = bRe;
    160.                         im = bIm;
    161.                         bRe = (re+bk[i])*zRe - im*zIm;
    162.                         bIm = (re+bk[i])*zIm + im*zRe;
    163.                 }
    164.                 bRe += bk[0];

    165.         // denominator
    166.                 aRe = 0;
    167.                 aIm = 0;
    168.                 for( int i=p; i>0; --i )
    169.                 {
    170.                         re = aRe;
    171.                         im = aIm;
    172.                         aRe = (re+ak[i])*zRe - im*zIm;
    173.                         aIm = (re+ak[i])*zIm + im*zRe;
    174.                 }
    175.                 aRe += ak[0];

    176.         // Power Spectrum Density
    177.                 numRe = aRe*bRe + aIm*bIm;
    178.                 numIm = aRe*bIm - aIm*bRe;
    179.                 den = aRe*aRe + aIm*aIm;
    180.                 Xre = numRe/(den);
    181.                 Xim = numIm/(den);
    182.         Xk[k] = sigma2 * (Xre*Xre + Xim*Xim);
    183.         }

    184.         return Xk;
    185. }
    复制代码


    测试代码:
    1. /*****************************************************************************
    2. *                            parametricpse_test.cpp
    3. *
    4. * parametric power spectrum estimation testing.
    5. *
    6. * Zhang Ming, 2010-11, Xi'an Jiaotong University.
    7. *****************************************************************************/


    8. #define BOUNDS_CHECK

    9. #include <iostream>
    10. #include <iomanip>
    11. #include <cstring>
    12. #include <random.h>
    13. #include <vectormath.h>
    14. #include <parametricpse.h>
    15. #include "engine.h"


    16. using namespace std;
    17. using namespace splab;


    18. typedef double  Type;
    19. const   int     N = 50;
    20. const   int     L = 200;
    21. const   int     yuleOrder = 4;
    22. const   int     burgOrder = 4;
    23. const   int     lplsOrder = 4;


    24. int main()
    25. {
    26.     /******************************* [ signal ] ******************************/
    27.     cout << setiosflags(ios::fixed) << setprecision(4);
    28.     int mfn = L/2+1;
    29.     Type amp1 = Type(1.0),
    30.          amp2 = Type(1.0);
    31.     Type f1 = Type(0.2),
    32.          f2 = Type(0.4);
    33.     Type sigma2, SNR;

    34.     Vector<Type> tn = linspace(Type(0), Type(N-1), N );
    35.     Vector<Type> sn = amp1*sin(TWOPI*f1*tn) + amp2*sin(TWOPI*f2*tn);
    36.     Vector<Type> wn = randn( 37, Type(0.0), Type(0.1), N );
    37.     Vector<Type> xn = sn + wn;
    38.     SNR = 20*log10(norm(sn)/norm(wn));
    39.     cout << "The SNR = " << SNR << endl << endl;

    40.     /********************************* [ PSD ] *******************************/
    41. //    Vector<Type> ak = yulewalkerPSE( xn, yuleOrder, sigma2 );
    42. //    cout << "The estimated AR coefficients by Youle-Walker method are: "
    43. //         << ak << endl;

    44. //    Vector<Type> ak = burgPSE( xn, burgOrder, sigma2 );
    45. //    cout << "The estimated AR coefficients by Burg method are: "
    46. //         << ak << endl;

    47.     Vector<Type> ak = fblplsPSE( xn, lplsOrder, sigma2 );
    48.     cout << "The estimated AR coefficients by Youle-Walker method are: "
    49.          << ak << endl;

    50.     cout << "The estimated variance is:   " << sigma2 << endl << endl;

    51.     Vector<Type> bk(1,Type(1.0));
    52.     Vector<Type> Px = armaPSD( ak, bk, sigma2, L );

    53.     /******************************** [ PLOT ] *******************************/
    54.         Engine *ep  = engOpen( NULL );
    55.         if( !ep )
    56.         {
    57.                 cerr << "Cannot open Matlab Engine!" << endl;
    58.                 exit(1);
    59.         }

    60.         mxArray *mxn = mxCreateDoubleMatrix( N, 1, mxREAL );
    61.     mxArray *mPx = mxCreateDoubleMatrix( mfn, 1, mxREAL );
    62.         memcpy( mxGetPr(mxn), xn, N*sizeof(Type) );
    63.         memcpy( mxGetPr(mPx), Px, mfn*sizeof(Type) );
    64.         engPutVariable( ep, "xn", mxn );
    65.         engPutVariable( ep, "Px", mPx );

    66.         const char *mCmd =  " figure('name','FBLPLS Method of Spectrum Estimation'); \
    67.         N = length(xn); mfn = length(Px); \
    68.         subplot(2,1,1); \
    69.             plot((0:N-1), xn); \
    70.             axis([0,N,min(xn),max(xn)]); \
    71.             title('(a)   Signal', 'FontSize',12); \
    72.             xlabel('Samples', 'FontSize',12); \
    73.             ylabel('Amplitude', 'FontSize',12); \
    74.         subplot(2,1,2); \
    75.             h = stem((0:mfn-1)/(mfn-1)/2, Px); \
    76.             axis([0,0.5,min(Px),max(Px)]); \
    77.             set(h,'MarkerFaceColor','blue'); \
    78.             set(gca, 'XTick', 0:0.1:0.5); \
    79.             grid on; \
    80.             title('(b)   Spectrum', 'FontSize',12); \
    81.             xlabel('Normalized Frequency ( f / fs )', 'FontSize',12); \
    82.             ylabel('Amplitude', 'FontSize',12); ";
    83.     engEvalString( ep, mCmd );

    84.         mxDestroyArray( mxn );
    85.         mxDestroyArray( mPx );
    86.     system( "pause" );
    87.         engClose(ep);

    88.         return 0;
    89. }
    复制代码


    运行结果:








    本文来自:http://my.oschina.net/zmjerry/blog/9480

    本帖子中包含更多资源

    您需要 登录 才可以下载或查看,没有帐号?立即注册

    x
    守望者AIR技术交流社区(www.airmyth.com)
    回复

    使用道具 举报

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

    本版积分规则

    
    关闭

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

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

    GMT+8, 2019-10-19 16:03 , Processed in 0.045794 second(s), 37 queries .

    守望者AIR

    守望者AIR技术交流社区

    本站成立于 2014年12月31日

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