- 积分
- 136080
- 注册时间
- 2014-12-27
- 最后登录
- 2024-3-8
- 在线时间
- 602 小时
- 威望
- 562
- 贡献
- 29
- 金币
- 52584
- 钢镚
- 1422
- 交易凭证
- 1
- 分享
- 0
- 精华
- 33
- 帖子
- 2094
- 主题
- 1742
TA的每日心情 | 擦汗 2018-4-10 15:18 |
---|
签到天数: 447 天 [LV.9]以坛为家II
超级版主
- 威望
- 562
- 贡献
- 29
- 金币
- 52584
- 钢镚
- 1422
|
头文件:
- /*
- * Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
- *
- * This program is free software; you can redistribute it and/or modify it
- * under the terms of the GNU General Public License as published by the
- * Free Software Foundation, either version 2 or any later version.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions are met:
- *
- * 1. Redistributions of source code must retain the above copyright notice,
- * this list of conditions and the following disclaimer.
- *
- * 2. Redistributions in binary form must reproduce the above copyright
- * notice, this list of conditions and the following disclaimer in the
- * documentation and/or other materials provided with the distribution.
- *
- * This program is distributed in the hope that it will be useful, but WITHOUT
- * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
- * more details. A copy of the GNU General Public License is available at:
- * http://www.fsf.org/licensing/licenses
- */
- /*****************************************************************************
- * toeplitz.h
- *
- * A Toeplitz matrix is defined by one row and one column. A symmetric
- * Toeplitz matrix is defined by just one row. Toeplitz generates Toeplitz
- * matrices given just the row or row and column description.
- *
- * Zhang Ming, 2010-11, Xi'an Jiaotong University.
- *****************************************************************************/
- #ifndef TOEPLITZ_H
- #define TOEPLITZ_H
- #include <vector.h>
- #include <matrix.h>
- namespace splab
- {
- template<typename Type> Matrix<Type> toeplitz( const Vector<Type>&,
- const Vector<Type>& );
- template<typename Type> Matrix<Type> toeplitz( const Vector<Type>& );
- #include <toeplitz-impl.h>
- }
- // namespace splab
- #endif
- // TOEPLITZ_H
复制代码- /*
- * Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
- *
- * This program is free software; you can redistribute it and/or modify it
- * under the terms of the GNU General Public License as published by the
- * Free Software Foundation, either version 2 or any later version.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions are met:
- *
- * 1. Redistributions of source code must retain the above copyright notice,
- * this list of conditions and the following disclaimer.
- *
- * 2. Redistributions in binary form must reproduce the above copyright
- * notice, this list of conditions and the following disclaimer in the
- * documentation and/or other materials provided with the distribution.
- *
- * This program is distributed in the hope that it will be useful, but WITHOUT
- * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
- * more details. A copy of the GNU General Public License is available at:
- * http://www.fsf.org/licensing/licenses
- */
- /*****************************************************************************
- * eigenanalysispse.h
- *
- * Eigenanalysis algorithms for spectrum estimation.
- *
- * The eigenanalysis algorithms perform eigen decomposition of the signal's
- * auto-correlation matrix (ACM) to estimate signal's frequency content. The
- * ACM can be decomposed into signal-subspace and noise-subspace, then use
- * the orthogonality of the tow subspaces to estimate the signal's spectrum.
- *
- * These algorithms, such as Pisarenko, MUSIC (MUltiple SIgnal Classification)
- * and ESPRIT (Estimation of Signal Parameters by Rotational Invariance
- * Techniques) are particularly suitable for signals that are the sum of
- * sinusoids with additive white Gaussian noise. The model order (the number
- * of complex exponential signal) is estimated by minimizing the MDL criterion.
- *
- * This file also provide the Capon's maximum likehood method or minimum
- * variance method for specturm estimation.
- *
- * Zhang Ming, 2010-11, Xi'an Jiaotong University.
- *****************************************************************************/
- #ifndef EIGENANALYSISPSE_H
- #define EIGENANALYSISPSE_H
- #include <toeplitz.h>
- #include <levinson.h>
- #include <linequs1.h>
- #include <svd.h>
- #include <evd.h>
- namespace splab
- {
- template<typename Type> Vector<Type> caponPSE( const Vector<Type>&,
- int, int );
- template<typename Type> Vector<Type> pisarenkoPSE( const Vector<Type>&,
- int, int, int );
- template<typename Type> Vector<Type> musicPSE( const Vector<Type>&,
- int, int, int );
- template<typename Type> Vector<Type> espritPSE( const Vector<Type>&,
- int, int );
- template<typename Type> int orderEst( const Vector<Type>&, int );
- #include <eigenanalysispse-impl.h>
- }
- // namespace splab
- #endif
- // EIGENANALYSISPSE_H
复制代码
实现文件:
- /*
- * Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
- *
- * This program is free software; you can redistribute it and/or modify it
- * under the terms of the GNU General Public License as published by the
- * Free Software Foundation, either version 2 or any later version.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions are met:
- *
- * 1. Redistributions of source code must retain the above copyright notice,
- * this list of conditions and the following disclaimer.
- *
- * 2. Redistributions in binary form must reproduce the above copyright
- * notice, this list of conditions and the following disclaimer in the
- * documentation and/or other materials provided with the distribution.
- *
- * This program is distributed in the hope that it will be useful, but WITHOUT
- * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
- * more details. A copy of the GNU General Public License is available at:
- * http://www.fsf.org/licensing/licenses
- */
- /*****************************************************************************
- * toeplitz-impl.h
- *
- * Implementationfor Toeplitz matrices generator.
- *
- * Zhang Ming, 2010-11, Xi'an Jiaotong University.
- *****************************************************************************/
- /**
- * Returns a nonsymmetric Toeplitz matrix T having cn as its first column and
- * rn as its first row. If the first elements of cn and rn are different, the
- * column element is used.
- */
- template <typename Type>
- Matrix<Type> toeplitz( const Vector<Type> &cn, const Vector<Type> &rn )
- {
- int M = cn.size(),
- N = rn.size();
- Matrix<Type> T(M,N);
- // main diagonal and below the main diagonal
- for( int d=0; d<M; ++d )
- for( int i=0; i<M-d; ++i )
- T[i+d][i] = cn[d];
- // above the main diagonal
- for( int d=1; d<N; ++d )
- for( int i=0; i<N-d; ++i )
- T[i][i+d] = rn[d];
- return T;
- }
- /**
- * Returns the symmetric or Hermitian Toeplitz matrix formed from vector rn,
- * where rn defines the first row of the matrix.
- */
- template <typename Type>
- Matrix<Type> toeplitz( const Vector<Type> &rn )
- {
- int N = rn.size();
- Matrix<Type> T(N,N);
- // main diagonal
- for( int i=0; i<N; ++i )
- T[i][i] = rn[0];
- // above and below the main diagonal
- for( int d=1; d<N; ++d )
- for( int i=0; i<N-d; ++i )
- T[i][i+d] = T[i+d][i] = rn[d];
- return T;
- }
复制代码- /*
- * Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
- *
- * This program is free software; you can redistribute it and/or modify it
- * under the terms of the GNU General Public License as published by the
- * Free Software Foundation, either version 2 or any later version.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions are met:
- *
- * 1. Redistributions of source code must retain the above copyright notice,
- * this list of conditions and the following disclaimer.
- *
- * 2. Redistributions in binary form must reproduce the above copyright
- * notice, this list of conditions and the following disclaimer in the
- * documentation and/or other materials provided with the distribution.
- *
- * This program is distributed in the hope that it will be useful, but WITHOUT
- * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
- * more details. A copy of the GNU General Public License is available at:
- * http://www.fsf.org/licensing/licenses
- */
- /*****************************************************************************
- * parametricpse-impl.h
- *
- * Implementation for eigenanalysis algorithms for spectrum estimation.
- *
- * Zhang Ming, 2010-11, Xi'an Jiaotong University.
- *****************************************************************************/
- /**
- * The Capon method (minimum variance method) for spectral estimation.
- * xn : input signal
- * M : the order of the covariance matrix
- * L : the number of estimated spectral samples
- * return : spectral estimates at L frequencies:
- * w = 0, 2*pi/L, ..., 2*pi(L-1)/L
- */
- template <typename Type>
- Vector<Type> caponPSE( const Vector<Type> &xn, int M, int L )
- {
- int N = xn.size();
- assert( M < N );
- Vector<Type> rn(M), Ere(M), Eim(M), Tre(M), Tim(M),
- tn(M), wn(M), Px(L);
- // auto-correlation matrix R
- for( int i=0; i<M; ++i )
- {
- tn[i] = Type(i);
- for( int k=0; k<N-i; ++k )
- rn[i] += xn[k+i]*xn[k];
- }
- rn /= Type(N-M);
- for( int k=0; k<L; ++k )
- {
- // real part and imaginary part of E(omega)
- wn = Type( TWOPI*k/L ) * tn;
- Ere = cos( wn );
- Eim = sin( wn );
- // Tre = inv(R)*Ere, Tim = inv(R)*Eim. Because R is Toeplitz, this can
- // be solved through Levinson algorithm.
- Tre = levinson( rn, Ere );
- Tim = levinson( rn, Eim );
- // denominator of spectrum
- Px[k] = dotProd(Ere,Tre) + dotProd(Eim,Tim);
- }
- return Type(M)/Px;
- }
- /**
- * The MUSIC method for spectral estimation.
- * xn : input signal
- * M : the order of the covariance matrix
- * p : the model order
- * L : the number of estimated spectral samples
- * return : spectral estimates (dB) at L frequencies:
- * w = 0, 2*pi/L, ..., 2*pi(L-1)/L
- */
- template <typename Type>
- Vector<Type> musicPSE( const Vector<Type> &xn, int M, int p, int L )
- {
- int N = xn.size();
- assert( M < N );
- assert( p < M );
- Type Tre, Tim, sum;
- Vector<Type> rm(M), tm(M),
- Ere(M), Eim(M), Vi(M),
- wk(M), Px(L);
- // auto-correlation matrix R
- for( int i=0; i<M; ++i )
- {
- tm[i] = Type(i);
- for( int k=0; k<N-i; ++k )
- rm[i] += xn[k+i]*xn[k];
- }
- rm /= Type(N-M);
- Matrix<Type> Rx = toeplitz(rm);
- SVD<Type> svd;
- svd.dec(Rx);
- Matrix<Type> U = svd.getU();
- for( int k=0; k<L; ++k )
- {
- // real part and imaginary part of E(omega)
- wk = Type( TWOPI*k/L ) * tm;
- Ere = cos( wk );
- Eim = sin( wk );
- // Tre = Ere*Vi, Tim = Eim*Vi;
- sum = 0;
- for( int i=p; i<M; ++i )
- {
- Vi = U.getColumn(i);
- Tre = dotProd( Ere, Vi );
- Tim = dotProd( Eim, Vi );
- sum += Tre*Tre + Tim*Tim;
- }
- // spectrum
- // Px[k] = 1 / ( sum );
- Px[k] = -10*log10(sum);
- }
- return Px;
- }
- /**
- * The Pisarenko method for spectral estimation.
- * xn : input signal
- * M : the order of the covariance matrix
- * p : the model order
- * L : the number of estimated spectral samples
- * return : spectral estimates (dB) at L frequencies:
- * w = 0, 2*pi/L, ..., 2*pi(L-1)/L
- */
- template <typename Type>
- Vector<Type> pisarenkoPSE( const Vector<Type> &xn, int M, int p, int L )
- {
- int N = xn.size();
- assert( M < N );
- assert( p < M );
- Type Tre, Tim;
- Vector<Type> rm(M), tm(M),
- Ere(M), Eim(M), Vi(M),
- wk(M), Px(L);
- // auto-correlation matrix R
- for( int i=0; i<M; ++i )
- {
- tm[i] = Type(i);
- for( int k=0; k<N-i; ++k )
- rm[i] += xn[k+i]*xn[k];
- }
- rm /= Type(N-M);
- Matrix<Type> Rx = toeplitz(rm);
- SVD<Type> svd;
- svd.dec(Rx);
- Matrix<Type> U = svd.getU();
- for( int k=0; k<L; ++k )
- {
- // real part and imaginary part of E(omega)
- wk = Type( TWOPI*k/L ) * tm;
- Ere = cos( wk );
- Eim = sin( wk );
- // Tre = Ere*Vi, Tim = Eim*Vi;
- Vi = U.getColumn(p);
- Tre = dotProd( Ere, Vi );
- Tim = dotProd( Eim, Vi );
- // spectrum
- // Px[k] = 1 / ( Tre*Tre + Tim*Tim );
- Px[k] = -10*log10(Tre*Tre+Tim*Tim);
- }
- return Px;
- }
- /**
- * The ESPRIT method for spectral estimation.
- * xn : input signal
- * M : the order of the covariance matrix
- * p : the model order
- * return : the esitmated frequency in the interval [-0.5,0.5]
- */
- template <typename Type>
- Vector<Type> espritPSE( Vector<Type> &xn, int M, int p )
- {
- int N = xn.size();
- assert( M < N );
- assert( p < M );
- Vector<Type> rm(M), fk(p);
- Matrix<Type> S1(M-1,p), S2(M-1,p);
- // get the auto-correlation matrix
- for( int i=0; i<M; ++i )
- for( int k=0; k<N-i; ++k )
- rm[i] += xn[k+i]*xn[k];
- rm /= Type(N-M);
- Matrix<Type> Rx = toeplitz(rm);
- // get the eigendecomposition of R, use svd because it sorts eigenvalues
- SVD<Type> svd;
- svd.dec(Rx);
- Matrix<Type> U = svd.getU();
- // compute S1 and S2
- for( int i=0; i<M-1; ++i )
- for( int j=0; j<p; ++j )
- {
- S1[i][j] = U[i][j];
- S2[i][j] = U[i+1][j];
- }
- // compute matrix Phi
- Matrix<Type> Phi = choleskySolver( trMult(S1,S1), trMult(S1,S2) );
- // compute eigenvalues of Phi
- EVD<Type> evd;
- evd.dec(Phi);
- Vector<Type> evRe = real( evd.getCD() );
- Vector<Type> evIM = imag( evd.getCD() );
- // compute normalized frequency in the interval [-0.5,0.5]
- for( int i=0; i<p; ++i )
- fk[i] = atan2( evIM[i], evRe[i] ) / Type(TWOPI);
- return fk;
- }
- /**
- * The model order estimation based on minimizing the MDL criterion.
- * xn : input signal
- * M : the order of the covariance matrix
- * return : p ---> the number of sinusoids
- */
- template <typename Type>
- int orderEst( const Vector<Type> &xn, int M )
- {
- int N = xn.size();
- assert( M < N );
- int p = 0;
- Type Gp, Ap, Ep, MDL, minMDL;
- Vector<Type> rm(M);
- // auto-correlation matrix R
- for( int i=0; i<M; ++i )
- for( int k=0; k<N-i; ++k )
- rm[i] += xn[k+i]*xn[k];
- rm /= Type(N-M);
- Matrix<Type> Rx = toeplitz(rm);
- SVD<Type> svd;
- svd.dec(Rx);
- Vector<Type> S = svd.getSV();
- minMDL = Type( 0.5*(M-1)*(M+1)*log10(1.0*N) );
- for( int i=0; i<M-2; ++i )
- {
- // compute MDL(i)
- Gp = Type(1);
- Ap = Type(0);
- Ep = Type( 0.5*i*(2*M-i)*log10(1.0*N) );
- for( int j=i+1; j<M; ++j )
- {
- Gp *= S[j];
- Ap += S[j];
- }
- Ap = pow( Ap/(M-i), Type(M-i) );
- MDL = -N*log10(Gp/Ap) + Ep;
- // find the minimum MDL(i)
- if( MDL < minMDL )
- {
- p = i;
- minMDL = MDL;
- }
- }
- return p;
- }
复制代码 测试代码:
- /*****************************************************************************
- * eigenanalysispse_test.cpp
- *
- * Eigenanalysis spectrum estimation testing.
- *
- * Zhang Ming, 2010-11, Xi'an Jiaotong University.
- *****************************************************************************/
- #define BOUNDS_CHECK
- #include <iostream>
- #include <iomanip>
- #include <cstring>
- #include <random.h>
- #include <vectormath.h>
- #include <eigenanalysispse.h>
- #include "engine.h"
- using namespace std;
- using namespace splab;
- typedef double Type;
- const int N = 200;
- const int M = 20;
- const int L = 200;
- int main()
- {
- /******************************* [ signal ] ******************************/
- cout << setiosflags(ios::fixed) << setprecision(4);
- int mfn = L/2+1;
- Type amp1 = Type(1.0),
- amp2 = Type(1.0);
- Type f1 = Type(0.2),
- f2 = Type(0.25);
- Type SNR;
- Vector<Type> tn = linspace(Type(0), Type(N-1), N );
- Vector<Type> sn = amp1*sin(Type(TWOPI)*f1*tn) + amp2*sin(Type(TWOPI)*f2*tn);
- Vector<Type> wn = randn( 37, Type(0.0), Type(0.2), N );
- Vector<Type> xn = sn + wn;
- SNR = 20*log10(norm(sn)/norm(wn));
- cout << "The SNR = " << SNR << endl << endl;
- /********************************* [ PSD ] *******************************/
- int p = orderEst( xn, M );
- p = p/2*2;
- cout << "The model order that minimizes the MDL criterion is: p = "
- << p << endl << endl;
- // Vector<Type> Px = caponPSE( xn, M, L );
- // Vector<Type> Px = pisarenkoPSE( xn, M, p, L );
- // Vector<Type> Px = musicPSE( xn, M, p, L );
- Vector<Type> fk = espritPSE( xn, M, p );
- Vector<Type> Px(mfn);
- for( int k=0; k<p; ++k )
- {
- int index = int( L*abs(fk[k]) + 0.5 );
- if( index != 0 && index != L/2 )
- Px[index] = Type(1.0);
- }
- /******************************** [ PLOT ] *******************************/
- Engine *ep = engOpen( NULL );
- if( !ep )
- {
- cerr << "Cannot open Matlab Engine!" << endl;
- exit(1);
- }
- mxArray *mxn = mxCreateDoubleMatrix( N, 1, mxREAL );
- mxArray *mPx = mxCreateDoubleMatrix( mfn, 1, mxREAL );
- memcpy( mxGetPr(mxn), xn, N*sizeof(Type) );
- memcpy( mxGetPr(mPx), Px, mfn*sizeof(Type) );
- engPutVariable( ep, "xn", mxn );
- engPutVariable( ep, "Px", mPx );
- const char *mCmd = " figure('name','FBLPLS Method of Spectrum Estimation'); \
- N = length(xn); mfn = length(Px); \
- subplot(2,1,1); \
- plot((0:N-1), xn); \
- axis([0,N,min(xn),max(xn)]); \
- title('(a) Signal', 'FontSize',12); \
- xlabel('Samples', 'FontSize',12); \
- ylabel('Amplitude', 'FontSize',12); \
- subplot(2,1,2); \
- h = stem((0:mfn-1)/(mfn-1)/2, Px); \
- axis([0,0.5,min(Px),max(Px)]); \
- set(h,'MarkerFaceColor','blue'); \
- set(gca, 'XTick', 0:0.05:0.5); \
- grid on; \
- title('(b) Spectrum', 'FontSize',12); \
- xlabel('Normalized Frequency ( f / fs )', 'FontSize',12); \
- ylabel('Amplitude', 'FontSize',12); ";
- engEvalString( ep, mCmd );
- mxDestroyArray( mxn );
- mxDestroyArray( mPx );
- system( "pause" );
- engClose(ep);
- return 0;
- }
复制代码 运行结果:
- The SNR = 14.3399
- The model order that minimizes the MDL criterion is: p = 4
- Process returned 0 (0x0) execution time : 0.062 s
- Press any key to continue.
复制代码 Capon估计法
MUSIC估计法
ESPRIT估计法
本文来自:http://my.oschina.net/zmjerry/blog/9584
|
本帖子中包含更多资源
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
|