经典谱估计算法(相关函数法,周期图就法,平滑周期图法)的C++实现

原创
2010/11/02 22:34
阅读数 5.8K

头文件:

/*
 * 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
 */


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


#ifndef CLASSICALPSE_H
#define CLASSICALPSE_H


#include <window.h>
#include <utilities.h>
#include <fft.h>
#include <correlation.h>


namespace splab
{

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

    template<typename Type> Vector<Type> periodogramPSE( const Vector<Type>&,
                                                         const Vector<Type>&,
                                                         int );
    template<typename Type> Vector<Type> bartlettPSE( const Vector<Type>&,
                                                      int, int );
    template<typename Type> Vector<Type> welchPSE( const Vector<Type>&,
                                                  const Vector<Type>&,
                                                  int, int );

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

    #include <classicalpse-impl.h>

}
// namespace splab


#endif
// CLASSICALSE_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
 */


/*****************************************************************************
 *                            classicalpse-impl.h
 *
 * Implementation for classical power spectrum estimatoin methods.
 *
 * Zhang Ming, 2010-11, Xi'an Jiaotong University.
 *****************************************************************************/


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


/**
 * The windowed periodogram power spectral estimator.
 * xn       : input signal
 * wn       : window function
 * L        : the number of power spectrum density samples
 * return   : spectral estimates at L frequencies:
 *            w = 0, 2*pi/L, ..., 2*pi(L-1)/L
 */
template <typename Type>
Vector<Type> periodogramPSE( const Vector<Type> &xn, const Vector<Type> &wn,
                             int L )
{
    int N = xn.size(),
	    M = wn.size();

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

    Vector<Type> wxn(L);
	if( L >= M )
	{
	    for( int i=0; i<M; ++i )
            wxn[i] = xn[i] * wn[i];
	}
	else
	{
        cerr << "The FFT points is smaller than the data points, ";
        cerr << "the data will be trucated to the FFT points!" << endl;
        for( int i=0; i<L; ++i )
            wxn[i] = xn[i] * wn[i];
	}

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


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


/**
 * The Welch method of power spectral estimation.
 * xn       : input signal
 * wn       : window function
 * K        : the number of subsequence
 * L        : the number of power spectrum density samples
 * return   : spectral estimates at L frequencies:
 *            w = 0, 2*pi/L, ..., 2*pi(L-1)/L
 */
template <typename Type>
Vector<Type> welchPSE( const Vector<Type> &xn, const Vector<Type> &wn,
                       int K, int L )
{
    int N = xn.size(),
	    M = wn.size();


	assert( M < N );
	assert( K < N );

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

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

	return phi/(S*P);
}


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

	assert( M <= N );

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

	Vector<Type> wrn(L);
	if( L >= M )
	{
	    for( int i=0; i<M; ++i )
            wrn[i] = Rxx[N-1+i] * wn[i];
	}
	else
	{
        cerr << "The FFT points is smaller than the data points, ";
        cerr << "the data will be trucated to the FFT points!" << endl;
        for( int i=0; i<L; ++i )
            wrn[i] = Rxx[N-1+i] * wn[i];
	}

	return Type(2)*real(fft(wrn)) - wrn[0];
}

测试代码:

/*****************************************************************************
 *                            classicalpse_test.cpp
 *
 * Classical power spectrum estimation testing.
 *
 * Zhang Ming, 2010-11, Xi'an Jiaotong University.
 *****************************************************************************/


#define BOUNDS_CHECK

#include <iostream>
#include <cstring>
#include <vectormath.h>
#include <classicalpse.h>
#include "engine.h"


using namespace std;
using namespace splab;


typedef double  Type;
const   int     N = 100;
const   int     K = 25;
const   int     M = 50;
const   int     L = 200;


int main()
{
    /******************************* [ signal ] ******************************/
    int mfn = L/2+1;
	Type amp1 = Type(1.0),
         amp2 = Type(0.5);
    Type f1 = Type(0.3),
         f2 = Type(0.4);
    Vector<Type> tn = linspace(Type(0), Type(N-1), N );
    Vector<Type> sn = amp1*sin(TWOPI*f1*tn) + amp2*sin(TWOPI*f2*tn);

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

	/********************************* [ PSD ] *******************************/
//	Vector<Type> Ps = correlogramPSE( sn, L );
//	Vector<Type> Ps = periodogramPSE( sn, wn, L );
//    Vector<Type> Ps = bartlettPSE( sn, M, L );
	Vector<Type> Ps = welchPSE( sn, wn, K, L );
//	Vector<Type> Ps = btPSE( sn, wn, L );

    /******************************** [ PLOT ] *******************************/
	Engine *ep  = engOpen( NULL );
	if( !ep )
	{
		cerr << "Cannot open Matlab Engine!" << endl;
		exit(1);
	}

	mxArray *msn = mxCreateDoubleMatrix( N, 1, mxREAL );
    mxArray *mPs = mxCreateDoubleMatrix( mfn, 1, mxREAL );
	memcpy( mxGetPr(msn), sn, N*sizeof(Type) );
	memcpy( mxGetPr(mPs), Ps, mfn*sizeof(Type) );
	engPutVariable( ep, "sn", msn );
	engPutVariable( ep, "Ps", mPs );

	const char *mCmd =  " figure('name','Welch Method of Spectrum Estimation'); \
        N = length(sn); mfn = length(Ps); \
        subplot(2,1,1); \
            plot((0:N-1), sn); \
            axis([0,N,min(sn),max(sn)]); \
            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, Ps); \
            axis([0,0.5,min(Ps),max(Ps)]); \
            set(h,'MarkerFaceColor','blue'); \
            set(gca, 'XTick', 0:0.1:0.5); \
            grid on; \
            title('(b)   Spectrum', 'FontSize',12); \
            xlabel('Normalized Frequency ( f / fs )', 'FontSize',12); \
            ylabel('Amplitude', 'FontSize',12); ";
    engEvalString( ep, mCmd );

	mxDestroyArray( msn );
	mxDestroyArray( mPs );
    system( "pause" );
	engClose(ep);

	return 0;
}

运行结果:

展开阅读全文
打赏
1
3 收藏
分享
加载中
打赏
3 评论
3 收藏
1
分享
返回顶部
顶部