守望者--AIR技术交流

 找回密码
 立即注册

QQ登录

只需一步,快速开始

扫一扫,访问微社区

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

[音频分析] 经典谱估计算法(相关函数法,周期图就法,平滑周期图...

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

    [LV.9]以坛为家II

    1742

    主题

    2094

    帖子

    13万

    积分

    超级版主

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

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

    开源英雄守望者

    发表于 2015-8-20 15:49:28 | 显示全部楼层 |阅读模式
    头文件:
    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. *                               classicalpse.h
    26. *
    27. * Classical Power Spectrum Estimation Mothods.
    28. *
    29. * The goal of spectral density estimation is to estimate the spectral density
    30. * of a random signal from a sequence of time samples of the signal. The
    31. * purpose of estimating the spectral density is to detect any periodicities
    32. * in the data, by observing peaks at the frequencies corresponding to these
    33. * periodicities.
    34. *
    35. * This file provides 5 usually used Classical-Specturm-Estimation methods:
    36. *    correlogram method,    periodogram method,
    37. *    smoothed periodogram method (Barteltt method and Welch method),
    38. *    and Blackman-Tukey method
    39. *
    40. * Zhang Ming, 2010-11, Xi'an Jiaotong University.
    41. *****************************************************************************/


    42. #ifndef CLASSICALPSE_H
    43. #define CLASSICALPSE_H


    44. #include <window.h>
    45. #include <utilities.h>
    46. #include <fft.h>
    47. #include <correlation.h>


    48. namespace splab
    49. {

    50.     template<typename Type> Vector<Type> correlogramPSE( const Vector<Type>&,
    51.                                                          int );

    52.     template<typename Type> Vector<Type> periodogramPSE( const Vector<Type>&,
    53.                                                          const Vector<Type>&,
    54.                                                          int );
    55.     template<typename Type> Vector<Type> bartlettPSE( const Vector<Type>&,
    56.                                                       int, int );
    57.     template<typename Type> Vector<Type> welchPSE( const Vector<Type>&,
    58.                                                   const Vector<Type>&,
    59.                                                   int, int );

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

    62.     #include <classicalpse-impl.h>

    63. }
    64. // namespace splab


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


    31. /**
    32. * The correlogram power spectral estimator.
    33. * xn       : input signal
    34. * L        : the number of power spectrum density samples
    35. * return   : spectral estimates at L frequencies:
    36. *            w = 0, 2*pi/L, ..., 2*pi(L-1)/L
    37. */
    38. template <typename Type>
    39. inline Vector<Type> correlogramPSE( const Vector<Type> &xn, int L )
    40. {
    41.     return periodogramPSE( xn, rectangle(xn.size(),Type(1)), L );
    42. }


    43. /**
    44. * The windowed periodogram power spectral estimator.
    45. * xn       : input signal
    46. * wn       : window function
    47. * L        : the number of power spectrum density samples
    48. * return   : spectral estimates at L frequencies:
    49. *            w = 0, 2*pi/L, ..., 2*pi(L-1)/L
    50. */
    51. template <typename Type>
    52. Vector<Type> periodogramPSE( const Vector<Type> &xn, const Vector<Type> &wn,
    53.                              int L )
    54. {
    55.     int N = xn.size(),
    56.             M = wn.size();

    57.         assert( M <= N );
    58.         if( M < N )
    59.         {
    60.             cerr << "The length of window is smaller than the length of data, ";
    61.         cerr << "the data will be trucated to the window length!" << endl;
    62.         }

    63.     Vector<Type> wxn(L);
    64.         if( L >= M )
    65.         {
    66.             for( int i=0; i<M; ++i )
    67.             wxn[i] = xn[i] * wn[i];
    68.         }
    69.         else
    70.         {
    71.         cerr << "The FFT points is smaller than the data points, ";
    72.         cerr << "the data will be trucated to the FFT points!" << endl;
    73.         for( int i=0; i<L; ++i )
    74.             wxn[i] = xn[i] * wn[i];
    75.         }

    76.         Vector<Type> absXk = abs( fft( wxn ) );
    77.         return absXk*absXk / Type(M);
    78. }


    79. /**
    80. * The Bartlett method of power spectral estimation.
    81. * xn       : input signal
    82. * M        : the length of subsequences
    83. * L        : the number of power spectrum density samples
    84. * return   : spectral estimates at L frequencies:
    85. *            w = 0, 2*pi/L, ..., 2*pi(L-1)/L
    86. */
    87. template <typename Type>
    88. inline Vector<Type> bartlettPSE( const Vector<Type> &xn, int M, int L )
    89. {
    90.     return welchPSE( xn, rectangle(M,Type(1)), M, L );
    91. }


    92. /**
    93. * The Welch method of power spectral estimation.
    94. * xn       : input signal
    95. * wn       : window function
    96. * K        : the number of subsequence
    97. * L        : the number of power spectrum density samples
    98. * return   : spectral estimates at L frequencies:
    99. *            w = 0, 2*pi/L, ..., 2*pi(L-1)/L
    100. */
    101. template <typename Type>
    102. Vector<Type> welchPSE( const Vector<Type> &xn, const Vector<Type> &wn,
    103.                        int K, int L )
    104. {
    105.     int N = xn.size(),
    106.             M = wn.size();


    107.         assert( M < N );
    108.         assert( K < N );

    109.     int S = ( N-M+K )/K;
    110.         Type P = sum( wn*wn ) / Type(M);

    111.         Vector<Type> phi(L);
    112.         for( int i=0; i<S; ++i )
    113.                 phi += periodogramPSE( wkeep(xn,M,i*K), wn, L );

    114.         return phi/(S*P);
    115. }


    116. /**
    117. * The Blackman-Tukey method of power spectral estimation.
    118. * The correlation function is obtained from the standard biased estimate.
    119. * xn       : input signal
    120. * wn       : window function
    121. * L        : the number of power spectrum density samples
    122. * return   : spectral estimates at L frequencies:
    123. *            w = 0, 2*pi/L, ..., 2*pi(L-1)/L
    124. */
    125. template <typename Type>
    126. Vector<Type> btPSE( const Vector<Type> &xn, const Vector<Type> &wn, int L )
    127. {
    128.     int N = xn.size(),
    129.             M = wn.size();

    130.         assert( M <= N );

    131.         Vector<Type> Rxx = fastCorr( xn, "biased" );

    132.         Vector<Type> wrn(L);
    133.         if( L >= M )
    134.         {
    135.             for( int i=0; i<M; ++i )
    136.             wrn[i] = Rxx[N-1+i] * wn[i];
    137.         }
    138.         else
    139.         {
    140.         cerr << "The FFT points is smaller than the data points, ";
    141.         cerr << "the data will be trucated to the FFT points!" << endl;
    142.         for( int i=0; i<L; ++i )
    143.             wrn[i] = Rxx[N-1+i] * wn[i];
    144.         }

    145.         return Type(2)*real(fft(wrn)) - wrn[0];
    146. }
    复制代码


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


    8. #define BOUNDS_CHECK

    9. #include <iostream>
    10. #include <cstring>
    11. #include <vectormath.h>
    12. #include <classicalpse.h>
    13. #include "engine.h"


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


    16. typedef double  Type;
    17. const   int     N = 100;
    18. const   int     K = 25;
    19. const   int     M = 50;
    20. const   int     L = 200;


    21. int main()
    22. {
    23.     /******************************* [ signal ] ******************************/
    24.     int mfn = L/2+1;
    25.         Type amp1 = Type(1.0),
    26.          amp2 = Type(0.5);
    27.     Type f1 = Type(0.3),
    28.          f2 = Type(0.4);
    29.     Vector<Type> tn = linspace(Type(0), Type(N-1), N );
    30.     Vector<Type> sn = amp1*sin(TWOPI*f1*tn) + amp2*sin(TWOPI*f2*tn);

    31.         /******************************** [ widow ] ******************************/
    32.         Vector<Type> wn = window( "Hamming", M, Type(1.0) );

    33.         /********************************* [ PSD ] *******************************/
    34. //        Vector<Type> Ps = correlogramPSE( sn, L );
    35. //        Vector<Type> Ps = periodogramPSE( sn, wn, L );
    36. //    Vector<Type> Ps = bartlettPSE( sn, M, L );
    37.         Vector<Type> Ps = welchPSE( sn, wn, K, L );
    38. //        Vector<Type> Ps = btPSE( sn, wn, L );

    39.     /******************************** [ PLOT ] *******************************/
    40.         Engine *ep  = engOpen( NULL );
    41.         if( !ep )
    42.         {
    43.                 cerr << "Cannot open Matlab Engine!" << endl;
    44.                 exit(1);
    45.         }

    46.         mxArray *msn = mxCreateDoubleMatrix( N, 1, mxREAL );
    47.     mxArray *mPs = mxCreateDoubleMatrix( mfn, 1, mxREAL );
    48.         memcpy( mxGetPr(msn), sn, N*sizeof(Type) );
    49.         memcpy( mxGetPr(mPs), Ps, mfn*sizeof(Type) );
    50.         engPutVariable( ep, "sn", msn );
    51.         engPutVariable( ep, "Ps", mPs );

    52.         const char *mCmd =  " figure('name','Welch Method of Spectrum Estimation'); \
    53.         N = length(sn); mfn = length(Ps); \
    54.         subplot(2,1,1); \
    55.             plot((0:N-1), sn); \
    56.             axis([0,N,min(sn),max(sn)]); \
    57.             title('(a)   Signal', 'FontSize',12); \
    58.             xlabel('Samples', 'FontSize',12); \
    59.             ylabel('Amplitude', 'FontSize',12); \
    60.         subplot(2,1,2); \
    61.             h = stem((0:mfn-1)/(mfn-1)/2, Ps); \
    62.             axis([0,0.5,min(Ps),max(Ps)]); \
    63.             set(h,'MarkerFaceColor','blue'); \
    64.             set(gca, 'XTick', 0:0.1:0.5); \
    65.             grid on; \
    66.             title('(b)   Spectrum', 'FontSize',12); \
    67.             xlabel('Normalized Frequency ( f / fs )', 'FontSize',12); \
    68.             ylabel('Amplitude', 'FontSize',12); ";
    69.     engEvalString( ep, mCmd );

    70.         mxDestroyArray( msn );
    71.         mxDestroyArray( mPs );
    72.     system( "pause" );
    73.         engClose(ep);

    74.         return 0;
    75. }
    复制代码


    运行结果:





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

    本帖子中包含更多资源

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

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

    使用道具 举报

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

    本版积分规则

    
    关闭

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

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

    GMT+8, 2019-8-19 01:02 , Processed in 0.044214 second(s), 37 queries .

    守望者AIR

    守望者AIR技术交流社区

    本站成立于 2014年12月31日

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