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

守望者--AIR技术交流

 找回密码
 立即注册

QQ登录

只需一步，快速开始

扫一扫，访问微社区

搜索
热搜: ANE FlasCC 炼金术
守望者--AIR技术交流»论坛 › 综合交流 › 专题资料 › 算法研究 › 经典谱估计算法（相关函数法，周期图就法，平滑周期图.. ...
返回列表 发新帖
查看: 2495|回复: 0
打印 上一主题 下一主题

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

[复制链接]
破晓
  • TA的每日心情
    擦汗
    2018-4-10 15:18

    • 签到天数: 447 天

    [LV.9]以坛为家II

    1742

    主题

    2094

    帖子

    13万

    积分

    超级版主

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

    威望
    562
    贡献
    29
    金币
    52695
    钢镚
    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. */


    26. /*****************************************************************************
    27. *                               classicalpse.h
    28. *
    29. * Classical Power Spectrum Estimation Mothods.
    30. *
    31. * The goal of spectral density estimation is to estimate the spectral density
    32. * of a random signal from a sequence of time samples of the signal. The
    33. * purpose of estimating the spectral density is to detect any periodicities
    34. * in the data, by observing peaks at the frequencies corresponding to these
    35. * periodicities.
    36. *
    37. * This file provides 5 usually used Classical-Specturm-Estimation methods:
    38. *    correlogram method,    periodogram method,
    39. *    smoothed periodogram method (Barteltt method and Welch method),
    40. *    and Blackman-Tukey method
    41. *
    42. * Zhang Ming, 2010-11, Xi'an Jiaotong University.
    43. *****************************************************************************/


    46. #ifndef CLASSICALPSE_H
    47. #define CLASSICALPSE_H


    50. #include <window.h>
    51. #include <utilities.h>
    52. #include <fft.h>
    53. #include <correlation.h>


    56. namespace splab
    57. {

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

    62.     template<typename Type> Vector<Type> periodogramPSE( const Vector<Type>&,
    63.                                                          const Vector<Type>&,
    64.                                                          int );
    65.     template<typename Type> Vector<Type> bartlettPSE( const Vector<Type>&,
    66.                                                       int, int );
    67.     template<typename Type> Vector<Type> welchPSE( const Vector<Type>&,
    68.                                                   const Vector<Type>&,
    69.                                                   int, int );

    71.     template<typename Type> Vector<Type> btPSE( const Vector<Type>&,
    72.                                                const Vector<Type>&, int );

    74.     #include <classicalpse-impl.h>

    76. }
    77. // namespace splab


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


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


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


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

    64.         assert( M <= N );
    65.         if( M < N )
    66.         {
    67.             cerr << "The length of window is smaller than the length of data, ";
    68.         cerr << "the data will be trucated to the window length!" << endl;
    69.         }

    71.     Vector<Type> wxn(L);
    72.         if( L >= M )
    73.         {
    74.             for( int i=0; i<M; ++i )
    75.             wxn[i] = xn[i] * wn[i];
    76.         }
    77.         else
    78.         {
    79.         cerr << "The FFT points is smaller than the data points, ";
    80.         cerr << "the data will be trucated to the FFT points!" << endl;
    81.         for( int i=0; i<L; ++i )
    82.             wxn[i] = xn[i] * wn[i];
    83.         }

    85.         Vector<Type> absXk = abs( fft( wxn ) );
    86.         return absXk*absXk / Type(M);
    87. }


    90. /**
    91. * The Bartlett method of power spectral estimation.
    92. * xn       : input signal
    93. * M        : the length of subsequences
    94. * L        : the number of power spectrum density samples
    95. * return   : spectral estimates at L frequencies:
    96. *            w = 0, 2*pi/L, ..., 2*pi(L-1)/L
    97. */
    98. template <typename Type>
    99. inline Vector<Type> bartlettPSE( const Vector<Type> &xn, int M, int L )
    100. {
    101.     return welchPSE( xn, rectangle(M,Type(1)), M, L );
    102. }


    105. /**
    106. * The Welch method of power spectral estimation.
    107. * xn       : input signal
    108. * wn       : window function
    109. * K        : the number of subsequence
    110. * L        : the number of power spectrum density samples
    111. * return   : spectral estimates at L frequencies:
    112. *            w = 0, 2*pi/L, ..., 2*pi(L-1)/L
    113. */
    114. template <typename Type>
    115. Vector<Type> welchPSE( const Vector<Type> &xn, const Vector<Type> &wn,
    116.                        int K, int L )
    117. {
    118.     int N = xn.size(),
    119.             M = wn.size();


    122.         assert( M < N );
    123.         assert( K < N );

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

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

    132.         return phi/(S*P);
    133. }


    136. /**
    137. * The Blackman-Tukey method of power spectral estimation.
    138. * The correlation function is obtained from the standard biased estimate.
    139. * xn       : input signal
    140. * wn       : window function
    141. * L        : the number of power spectrum density samples
    142. * return   : spectral estimates at L frequencies:
    143. *            w = 0, 2*pi/L, ..., 2*pi(L-1)/L
    144. */
    145. template <typename Type>
    146. Vector<Type> btPSE( const Vector<Type> &xn, const Vector<Type> &wn, int L )
    147. {
    148.     int N = xn.size(),
    149.             M = wn.size();

    151.         assert( M <= N );

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

    155.         Vector<Type> wrn(L);
    156.         if( L >= M )
    157.         {
    158.             for( int i=0; i<M; ++i )
    159.             wrn[i] = Rxx[N-1+i] * wn[i];
    160.         }
    161.         else
    162.         {
    163.         cerr << "The FFT points is smaller than the data points, ";
    164.         cerr << "the data will be trucated to the FFT points!" << endl;
    165.         for( int i=0; i<L; ++i )
    166.             wrn[i] = Rxx[N-1+i] * wn[i];
    167.         }

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


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


    10. #define BOUNDS_CHECK

    12. #include <iostream>
    13. #include <cstring>
    14. #include <vectormath.h>
    15. #include <classicalpse.h>
    16. #include "engine.h"


    19. using namespace std;
    20. using namespace splab;


    23. typedef double  Type;
    24. const   int     N = 100;
    25. const   int     K = 25;
    26. const   int     M = 50;
    27. const   int     L = 200;


    30. int main()
    31. {
    32.     /******************************* [ signal ] ******************************/
    33.     int mfn = L/2+1;
    34.         Type amp1 = Type(1.0),
    35.          amp2 = Type(0.5);
    36.     Type f1 = Type(0.3),
    37.          f2 = Type(0.4);
    38.     Vector<Type> tn = linspace(Type(0), Type(N-1), N );
    39.     Vector<Type> sn = amp1*sin(TWOPI*f1*tn) + amp2*sin(TWOPI*f2*tn);

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

    44.         /********************************* [ PSD ] *******************************/
    45. //        Vector<Type> Ps = correlogramPSE( sn, L );
    46. //        Vector<Type> Ps = periodogramPSE( sn, wn, L );
    47. //    Vector<Type> Ps = bartlettPSE( sn, M, L );
    48.         Vector<Type> Ps = welchPSE( sn, wn, K, L );
    49. //        Vector<Type> Ps = btPSE( sn, wn, L );

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

    59.         mxArray *msn = mxCreateDoubleMatrix( N, 1, mxREAL );
    60.     mxArray *mPs = mxCreateDoubleMatrix( mfn, 1, mxREAL );
    61.         memcpy( mxGetPr(msn), sn, N*sizeof(Type) );
    62.         memcpy( mxGetPr(mPs), Ps, mfn*sizeof(Type) );
    63.         engPutVariable( ep, "sn", msn );
    64.         engPutVariable( ep, "Ps", mPs );

    66.         const char *mCmd =  " figure('name','Welch Method of Spectrum Estimation'); \
    67.         N = length(sn); mfn = length(Ps); \
    68.         subplot(2,1,1); \
    69.             plot((0:N-1), sn); \
    70.             axis([0,N,min(sn),max(sn)]); \
    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, Ps); \
    76.             axis([0,0.5,min(Ps),max(Ps)]); \
    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 );

    85.         mxDestroyArray( msn );
    86.         mxDestroyArray( mPs );
    87.     system( "pause" );
    88.         engClose(ep);

    90.         return 0;
    91. }
    复制代码


    运行结果：





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

    本帖子中包含更多资源

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

    x
    分享到:  QQ好友和群QQ好友和群 QQ空间QQ空间 腾讯微博腾讯微博 腾讯朋友腾讯朋友 微信微信
    收藏收藏 分享分享 支持支持 反对反对 微信

    相关帖子

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

    使用道具 举报

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

    本版积分规则

    
    关闭

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

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

    GMT+8, 2024-4-20 14:05 , Processed in 0.054156 second(s), 36 queries .

    守望者AIR

    守望者AIR技术交流社区

    本站成立于 2014年12月31日

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