头文件：
/*
* 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;
}运行结果：
677
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
