头文件：
/*
* 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估计法
572
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
