设为首页收藏本站
开启辅助访问 打卡签到论坛帮助切换到窄版

守望者--AIR技术交流

 找回密码
 立即注册

QQ登录

只需一步，快速开始

扫一扫，访问微社区

搜索
热搜: ANE FlasCC 炼金术
守望者--AIR技术交流»论坛 › 综合交流 › 专题资料 › 算法研究 › 基于AR模型谱估计算法（Yule-Walker方法与Burg方法）的C ...
返回列表 发新帖
查看: 1823|回复: 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
    金币
    52683
    钢镚
    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. */


    26. /*****************************************************************************
    27. *                                  levinson.h
    28. *
    29. * If the coefficient matrix of a linear equations is Toeplitz, then it can
    30. * be solved in a high computational efficiency way through Levinson-Durbin
    31. * algorithm. The subroutiones in this file will be used for solving Wiener
    32. * -Hopf equeations in Wiener filtring and Yule- Walker equations in
    33. * parametric spectrum estimation, respectively.
    34. *
    35. * Zhang Ming, 2010-11, Xi'an Jiaotong University.
    36. *****************************************************************************/


    39. #ifndef LEVINSON_H
    40. #define LEVINSON_H


    43. #include <vector.h>


    46. namespace splab
    47. {

    49.     template<typename Type>
    50.     Vector<Type> levinson( const Vector<Type>&, const Vector<Type>& );

    52.     template<typename Type>
    53.     Vector<Type> levinson( const Vector<Type>&, Type& );


    56.     #include <levinson-impl.h>

    58. }
    59. // namespace splab


    62. #endif
    63. // 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. */


    26. /*****************************************************************************
    27. *                               parametricpse.h
    28. *
    29. * Parametric Power Spectrum Estimation Mothods.
    30. *
    31. * Techniques for spectrum estimation can generally be divided into parametric
    32. * (such as classical spectrum estimation) and non-parametric methods. The
    33. * parametric approaches assume that the underlying stationary stochastic
    34. * process has a certain structure which can be described using a small number
    35. * of parameters (for example, auto-regressive model). In these approaches,
    36. * the task is to estimate the parameters of the model that describes the
    37. * stochastic process.
    38. *
    39. * The widely used model is AR model, so this file provides three subroutines
    40. * to estimate the parameter of AR model, they are Yule-Walker method, Burg's
    41. * recursive mothod and forward-and-backward linear prediction least square
    42. * method.
    43. *
    44. * Zhang Ming, 2010-11, Xi'an Jiaotong University.
    45. *****************************************************************************/


    48. #ifndef PARAMETRICPSE_H
    49. #define PARAMETRICPSE_H


    52. #include <vector.h>
    53. #include <fft.h>
    54. #include <correlation.h>
    55. #include <levinson.h>
    56. #include <linequs1.h>


    59. namespace splab
    60. {

    62.     template<typename Type> Vector<Type> yulewalkerPSE( const Vector<Type>&,
    63.                                                         int, Type& );
    64.     template<typename Type> Vector<Type> burgPSE( const Vector<Type>&,
    65.                                                   int, Type& );
    66.     template<typename Type> Vector<Type> fblplsPSE( const Vector<Type>&,
    67.                                                     int, Type& );

    69.     template<typename Type> Vector<Type> armaPSD( const Vector<Type>&,
    70.                                                   const Vector<Type>&,
    71.                                                   const Type&, int );


    74.     #include <parametricpse-impl.h>

    76. }
    77. // namespace splab


    80. #endif
    81. // 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. */


    26. /*****************************************************************************
    27. *                               levinson-impl.h
    28. *
    29. * Implementationfor Levinson-Durbin alogrithm.
    30. *
    31. * Zhang Ming, 2010-11, Xi'an Jiaotong University.
    32. *****************************************************************************/


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

    45.     int n = t.size();
    46.         Type alpha, beta, q, c, omega;
    47.         Vector<Type> y(n), yy(n), x(n);

    49.         alpha = t[0];
    50.         if( abs(alpha) < EPS )
    51.         {
    52.                 cerr << "The matrix is ill-conditioned!" << endl;
    53.                 return x;
    54.         }
    55.         y[0] = 1;
    56.         x[0] = b[0] / alpha;

    58.         for( int k=1; k<n; ++k )
    59.         {
    60.                 q = 0;
    61.                 beta = 0;
    62.                 for( int j=0; j<k; ++j )
    63.                 {
    64.                         q += x[j] * t[k-j];
    65.                         beta += y[j] * t[j+1];
    66.                 }
    67.                 c = -beta / alpha;

    69.                 yy[0] = c * y[k-1];
    70.                 y[k] = y[k-1];
    71.                 for( int i=1; i<k; ++i )
    72.                         yy[i] = y[i-1] + c*y[k-i-1];
    73.                 yy[k] = y[k-1];

    75.                 alpha += c*beta;
    76.                 if( abs(alpha) < EPS )
    77.                 {
    78.                         cerr << "The matrix is ill-conditioned!" << endl;
    79.                         return x;
    80.                 }

    82.                 omega = (b[k]-q) / alpha;
    83.                 for( int i=0; i<k; ++i )
    84.                 {
    85.                         x[i] += omega*yy[i];
    86.                         y[i] = yy[i];
    87.                 }
    88.                 x[k] = omega*y[k];
    89.         }

    91.         return x;
    92. }


    95. /**
    96. * Levinson-Durbin algorithm for solving Youle-Walker equations.
    97. * rn       : r(0), r(1), ..., r(p)
    98. * sigma2   : the variance of exciting white noise
    99. */
    100. template <typename Type>
    101. Vector<Type> levinson( const Vector<Type> &rn, Type &sigma2 )
    102. {
    103.     int p = rn.size()-1;
    104.     Type tmp;
    105.     Vector<Type> ak(p+1), akPrev(p+1);

    107.     ak[0] = Type(1.0);
    108.     sigma2 = rn[0];
    109.     ak[1] = -rn[1]/sigma2;
    110.     sigma2 *= 1 - ak[1]*ak[1];

    112.     for( int k=2; k<=p; ++k )
    113.     {
    114.         tmp = 0;
    115.         for( int i=0; i<k; ++i )
    116.             tmp += ak[i]*rn[k-i];
    117.         ak[k] = -tmp/sigma2;

    119.         for( int i=1; i<k; ++i )
    120.             akPrev[i] = ak[i] + ak[k]*ak[k-i];
    121.         for( int i=1; i<k; ++i )
    122.             ak[i] = akPrev[i];

    124.         sigma2 *= 1 - ak[k]*ak[k];
    125.     }

    127.         return ak;
    128. }
    复制代码
    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. */


    26. /*****************************************************************************
    27. *                            parametricpse-impl.h
    28. *
    29. * Implementation for parametric power spectrum estimatoin methods.
    30. *
    31. * Zhang Ming, 2010-11, Xi'an Jiaotong University.
    32. *****************************************************************************/


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

    47.     assert( p <= N );

    49.     Vector<Type> rn(p+1);
    50.     for( int i=0; i<=p; ++i )
    51.         for( int k=0; k<N-i; ++k )
    52.             rn[i] += xn[k+i]*xn[k];
    53.     rn /= Type(N);

    55.     return levinson( rn, sigma2 );
    56. }


    59. /**
    60. * The Burg method for AR power spectral estimation.
    61. * xn       : input signal
    62. * p        : the AR model order
    63. * sigma2   : the variance of exciting white noise
    64. * return   : coefficients of AR model --- a(0), a(1), ..., a(p)
    65. */
    66. template <typename Type>
    67. Vector<Type> burgPSE( const Vector<Type> &xn, int p, Type &sigma2 )
    68. {
    69.     int N = xn.size();
    70.     Type numerator, denominator;
    71.     Vector<Type> ak(p+1), akPrev(p+1), ef(N), eb(N);

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

    76.     for( int i=1; i<N; ++i )
    77.     {
    78.         ef[i] = xn[i];
    79.         eb[i-1] = xn[i-1];
    80.     }

    82.     for( int k=1; k<=p; ++k )
    83.     {
    84.         numerator = 0;
    85.         denominator = 0;
    86.         for( int i=k; i<N; ++i )
    87.         {
    88.             numerator += ef[i]*eb[i-1];
    89.             denominator += ef[i]*ef[i] + eb[i-1]*eb[i-1];
    90.         }
    91.         ak[k] = -2*numerator/denominator;

    93.         for( int i=1; i<k; ++i )
    94.             akPrev[i] = ak[i] + ak[k]*ak[k-i];
    95.         for( int i=1; i<k; ++i )
    96.             ak[i] = akPrev[i];

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

    100.         for( int i=N-1; i>k; --i )
    101.         {
    102.             ef[i] = ef[i] + ak[k]*eb[i-1];
    103.             eb[i-1] = eb[i-2] + ak[k]*ef[i-1];
    104.         }
    105.     }

    107.     return ak;
    108. }


    111. /**
    112. * Forward and backward linear prediction least square method for
    113. * AR power spectral estimation.
    114. * xn       : input signal
    115. * p        : the AR model order
    116. * sigma2   : the variance of exciting white noise
    117. * return   : coefficients of AR model --- a(0), a(1), ..., a(p)
    118. */
    119. template <typename Type>
    120. Vector<Type> fblplsPSE( const Vector<Type> &xn, int p, Type &sigma2 )
    121. {
    122.     int N = xn.size(),
    123.         M = 2*(N-p);

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

    128.     Matrix<Type> X(M,p+1);
    129. //    for( int i=0; i<=p; ++i )
    130. //    {
    131. //        for( int j=0; j<N-p; ++j )
    132. //            X[j][i] = xn[i+j];
    133. //        for( int j=p; j<N; ++j )
    134. //            X[j][i] = xn[j-i];
    135. //    }
    136.     for( int i=0; i<N-p; ++i )
    137.         for( int j=0; j<=p; ++j )
    138.             X[i][j] = xn[i+j];
    139.     for( int i=p; i<N; ++i )
    140.         for( int j=0; j<=p; ++j )
    141.             X[i][j] = xn[i-j];

    143.     Matrix<Type> Rp = trMult(X,X) / Type(M);
    144.     Vector<Type> ak = luSolver( Rp, u );
    145.     sigma2 = 1/ak[0];
    146.     ak *= sigma2;

    148.     return ak;
    149. }


    152. /**
    153. * The Bartlett method of power spectral estimation.
    154. * ak       : AR coefficients
    155. * bk       : MA coefficients
    156. * sigma2   : the variance of exciting white noise
    157. * L        : the points number of PSD
    158. * return   : spectral density at L frequencies:
    159. *            w = 0, 2*pi/L, ..., 2*pi(L-1)/L
    160. */
    161. template <typename Type>
    162. Vector<Type> armaPSD( const Vector<Type> &ak, const Vector<Type> &bk,
    163.                       const Type &sigma2, int L )
    164. {
    165.     int p = ak.size()-1,
    166.         q = bk.size()-1;
    167.     Vector<Type> Xk(L);

    169.     Type zRe, zIm, aRe, aIm, bRe, bIm,
    170.          Xre, Xim,
    171.          re, im;
    172.         Type omega,
    173.          den, numRe, numIm;

    175.         for( int k=0; k<L; ++k )
    176.         {
    177.                 omega = Type(TWOPI*k/L);
    178.                 zRe = cos(-omega);
    179.                 zIm = sin(-omega);

    181.         // numerator
    182.                 bRe = 0;
    183.                 bIm = 0;
    184.                 for( int i=q; i>0; --i )
    185.                 {
    186.                         re = bRe;
    187.                         im = bIm;
    188.                         bRe = (re+bk[i])*zRe - im*zIm;
    189.                         bIm = (re+bk[i])*zIm + im*zRe;
    190.                 }
    191.                 bRe += bk[0];

    193.         // denominator
    194.                 aRe = 0;
    195.                 aIm = 0;
    196.                 for( int i=p; i>0; --i )
    197.                 {
    198.                         re = aRe;
    199.                         im = aIm;
    200.                         aRe = (re+ak[i])*zRe - im*zIm;
    201.                         aIm = (re+ak[i])*zIm + im*zRe;
    202.                 }
    203.                 aRe += ak[0];

    205.         // Power Spectrum Density
    206.                 numRe = aRe*bRe + aIm*bIm;
    207.                 numIm = aRe*bIm - aIm*bRe;
    208.                 den = aRe*aRe + aIm*aIm;
    209.                 Xre = numRe/(den);
    210.                 Xim = numIm/(den);
    211.         Xk[k] = sigma2 * (Xre*Xre + Xim*Xim);
    212.         }

    214.         return Xk;
    215. }
    复制代码


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


    10. #define BOUNDS_CHECK

    12. #include <iostream>
    13. #include <iomanip>
    14. #include <cstring>
    15. #include <random.h>
    16. #include <vectormath.h>
    17. #include <parametricpse.h>
    18. #include "engine.h"


    21. using namespace std;
    22. using namespace splab;


    25. typedef double  Type;
    26. const   int     N = 50;
    27. const   int     L = 200;
    28. const   int     yuleOrder = 4;
    29. const   int     burgOrder = 4;
    30. const   int     lplsOrder = 4;


    33. int main()
    34. {
    35.     /******************************* [ signal ] ******************************/
    36.     cout << setiosflags(ios::fixed) << setprecision(4);
    37.     int mfn = L/2+1;
    38.     Type amp1 = Type(1.0),
    39.          amp2 = Type(1.0);
    40.     Type f1 = Type(0.2),
    41.          f2 = Type(0.4);
    42.     Type sigma2, SNR;

    44.     Vector<Type> tn = linspace(Type(0), Type(N-1), N );
    45.     Vector<Type> sn = amp1*sin(TWOPI*f1*tn) + amp2*sin(TWOPI*f2*tn);
    46.     Vector<Type> wn = randn( 37, Type(0.0), Type(0.1), N );
    47.     Vector<Type> xn = sn + wn;
    48.     SNR = 20*log10(norm(sn)/norm(wn));
    49.     cout << "The SNR = " << SNR << endl << endl;

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

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

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

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

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

    69.     /******************************** [ PLOT ] *******************************/
    70.         Engine *ep  = engOpen( NULL );
    71.         if( !ep )
    72.         {
    73.                 cerr << "Cannot open Matlab Engine!" << endl;
    74.                 exit(1);
    75.         }

    77.         mxArray *mxn = mxCreateDoubleMatrix( N, 1, mxREAL );
    78.     mxArray *mPx = mxCreateDoubleMatrix( mfn, 1, mxREAL );
    79.         memcpy( mxGetPr(mxn), xn, N*sizeof(Type) );
    80.         memcpy( mxGetPr(mPx), Px, mfn*sizeof(Type) );
    81.         engPutVariable( ep, "xn", mxn );
    82.         engPutVariable( ep, "Px", mPx );

    84.         const char *mCmd =  " figure('name','FBLPLS Method of Spectrum Estimation'); \
    85.         N = length(xn); mfn = length(Px); \
    86.         subplot(2,1,1); \
    87.             plot((0:N-1), xn); \
    88.             axis([0,N,min(xn),max(xn)]); \
    89.             title('(a)   Signal', 'FontSize',12); \
    90.             xlabel('Samples', 'FontSize',12); \
    91.             ylabel('Amplitude', 'FontSize',12); \
    92.         subplot(2,1,2); \
    93.             h = stem((0:mfn-1)/(mfn-1)/2, Px); \
    94.             axis([0,0.5,min(Px),max(Px)]); \
    95.             set(h,'MarkerFaceColor','blue'); \
    96.             set(gca, 'XTick', 0:0.1:0.5); \
    97.             grid on; \
    98.             title('(b)   Spectrum', 'FontSize',12); \
    99.             xlabel('Normalized Frequency ( f / fs )', 'FontSize',12); \
    100.             ylabel('Amplitude', 'FontSize',12); ";
    101.     engEvalString( ep, mCmd );

    103.         mxDestroyArray( mxn );
    104.         mxDestroyArray( mPx );
    105.     system( "pause" );
    106.         engClose(ep);

    108.         return 0;
    109. }
    复制代码


    运行结果：








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

    本帖子中包含更多资源

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

    x

    相关帖子

    守望者AIR技术交流社区（www.airmyth.com）
    回复

    使用道具 举报

    返回列表 发新帖
    高级模式
    B Color Image Link Quote Code Smilies
    您需要登录后才可以回帖 登录 | 立即注册

    本版积分规则

    
    关闭

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

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

    GMT+8, 2024-4-17 01:39 , Processed in 0.054059 second(s), 37 queries .

    守望者AIR

    守望者AIR技术交流社区

    本站成立于 2014年12月31日

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