相关及其快速算法的C++实现

2019独角兽企业重金招聘Python工程师标准>>>

头文件：

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


/*****************************************************************************
 *                                 correlation.h
 *
 * The routines in this file estimate the cross-correlation sequence of a
 * random process. Autocorrelation is handled as a special case.
 *
 * c = corr(x,y,opt) returns the cross-correlation sequence in a length
 * 2*N-1 vector, where x and y are length N vectors (N>1). If x and y are
 * not the same length, the shorter vector is zero-padded to the length of
 * the longer vector.
 *
 * The parameter "opt" specifies a normalization option for the cross-
 * correlation, where 'opt' is
 * "none"       : to use the raw, unscaled cross-correlations (default);
 * "biased"     : biased estimate of the cross-correlation function;
 * "unbiased"   : unbiased estimate of the cross-correlation function.
 *
 * We use FFT computing the auto-corelation and cross-corelation functions
 * based on fallowing facts: for real functions,
 * R1[x(t),y(t)] = sum{ x(u)*y(u-t) } = Conv[x(t),y(-t)]
 * R2[x(t),y(t)] = sum{ x(u)*y(u+t) } = Conv[x(-t),y(t)]
 * And here we use the first defination.
 *
 * Zhang Ming, 2010-10, Xi'an Jiaotong University.
 *****************************************************************************/


#ifndef CORRELATION_H
#define CORRELATION_H


#include <convolution.h>
#include <utilities.h>


namespace splab
{

    template<typename Type> Vector<Type> corr( const Vector<Type>&,
                                               const string &opt="none" );
    template<typename Type> Vector<Type> corr( const Vector<Type>&,
                                               const Vector<Type>&,
                                               const string &opt="none" );

    template<typename Type> Vector<Type> fastCorr( const Vector<Type>&,
                                                   const string &opt="none" );
    template<typename Type> Vector<Type> fastCorr( const Vector<Type>&,
                                                   const Vector<Type>&,
                                                   const string &opt="none" );

    template<typename Type> static void biasedProcessing( Vector<Type> &,
                                                          const string &opt );


    #include <correlation-impl.h>

}
// namespace splab


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


/*****************************************************************************
 *                             correlation-impl.h
 *
 * Implementation for linear correlation.
 *
 * Zhang Ming, 2010-10, Xi'an Jiaotong University.
 *****************************************************************************/


/**
 * Auto-correlation by defination in time domain.
 */
template<typename Type>
inline Vector<Type> corr( const Vector<Type> &xn, const string &opt )
{
    Vector<Type> rn = conv( xn, reverse(xn) );

    biasedProcessing( rn, opt );

    return rn;
}


/**
 * Cross-correlation by defination in time domain.
 */
template<typename Type>
inline Vector<Type> corr( const Vector<Type> &xn, const Vector<Type> &yn,
                          const string &opt )
{
    int N = xn.size(),
        d = N - yn.size();
    Vector<Type> rn;

    if( d > 0 )
        rn = conv( xn, reverse(wextend(yn,d,"right","zpd")) );
    else if( d < 0 )
    {
        N -= d;
        rn = conv( wextend(xn,-d,"right","zpd"), reverse(yn) );
    }
    else
        rn = conv( xn, reverse(yn) );

    biasedProcessing( rn, opt );

    return rn;
}


/**
 * Fast auto-correlation by using FFT.
 */
template<typename Type>
inline Vector<Type> fastCorr( const Vector<Type> &xn, const string &opt )
{
    Vector<Type> rn = fastConv( xn, reverse(xn) );

    biasedProcessing( rn, opt );

    return rn;
}


/**
 * Fast cross-correlation by using FFT.
 */
template<typename Type>
inline Vector<Type> fastCorr( const Vector<Type> &xn, const Vector<Type> &yn,
                              const string &opt )
{
    int N = xn.size(),
        d = N - yn.size();
    Vector<Type> rn;

    if( d > 0 )
        rn = fastConv( xn, reverse(wextend(yn,d,"right","zpd")) );
    else if( d < 0 )
    {
        N -= d;
        rn = fastConv( wextend(xn,-d,"right","zpd"), reverse(yn) );
    }
    else
        rn = fastConv( xn, reverse(yn) );

    biasedProcessing( rn, opt );

    return rn;
}


/**
 * Biase processing for correlation.
 */
template<typename Type>
static void biasedProcessing( Vector<Type> &rn, const string &opt )
{
    int N = (rn.size()+1) / 2;

    if( opt == "biased" )
        rn /= Type(N);
    else if( opt == "unbiased" )
    {
        int mid = N-1;
        rn[mid] /= N;
        for( int i=1; i<N; ++i )
        {
            rn[mid+i] /= (N-i);
            rn[mid-i] /= (N-i);
        }
    }
}

测试代码：

/*****************************************************************************
 *                             correlation_test.cpp
 *
 * Correlation testing.
 *
 * Zhang Ming, 2010-10, Xi'an Jiaotong University.
 *****************************************************************************/


#define BOUNDS_CHECK

#include <iostream>
#include <iomanip>
#include <correlation.h>


using namespace std;
using namespace splab;


typedef double  Type;
const   int     M = 3;
const   int     N = 5;


int main()
{
    Vector<Type> xn( M ), yn( N );

    for( int i=0; i<M; ++i )
        xn[i] = i;
    for( int i=0; i<N; ++i )
        yn[i] = i-N/2;

    cout << setiosflags(ios::fixed) << setprecision(4);
    cout << "xn:   " << xn << endl;
    cout << "yn:   " << yn << endl;

    // auto and cross correlation functions
    cout << "auto-correlation of xn:   " << corr(xn) << endl;
    cout << "biased auto-correlation of xn:   " << corr(xn,"biased") << endl;
    cout << "unbiased auto-correlation of xn:   " << corr(xn,"unbiased") << endl;

    cout << "cross-correlation of xn and yn:   " << corr(xn,yn) << endl;
    cout << "cross-correlation of yn and xn:   " << corr(yn,xn) << endl;

    cout << "biased cross-correlation of xn and yn:   "
         << corr(xn,yn,"biased") << endl;
    cout << "biased cross-correlation of yn and xn:   "
         << corr(yn,xn,"biased") << endl;

    cout << "unbiased cross-correlation of xn and yn:   "
         << corr(xn,yn,"unbiased") << endl;
    cout << "unbiased cross-correlation of yn and xn:   "
         << corr(yn,xn,"unbiased") << endl;

    // fast auto and cross correlation functions
    cout << "fast auto-correlation of xn:   " << fastCorr(xn) << endl;
    cout << "fast biased auto-correlation of xn:   " << fastCorr(xn,"biased") << endl;
    cout << "fast unbiased auto-correlation of xn:   " << fastCorr(xn,"unbiased") << endl;

    cout << "fast cross-correlation of xn and yn:   " << fastCorr(xn,yn) << endl;
    cout << "fast cross-correlation of yn and xn:   " << fastCorr(yn,xn) << endl;

    cout << "fast biased cross-correlation of xn and yn:   "
         << fastCorr(xn,yn,"biased") << endl;
    cout << "fast biased cross-correlation of yn and xn:   "
         << fastCorr(yn,xn,"biased") << endl;

    cout << "fast unbiased cross-correlation of xn and yn:   "
         << fastCorr(xn,yn,"unbiased") << endl;
    cout << "fast unbiased cross-correlation of yn and xn:   "
         << fastCorr(yn,xn,"unbiased") << endl;

    return 0;
}

运行结果：

xn:   size: 3 by 1
0.0000
1.0000
2.0000

yn:   size: 5 by 1
-2.0000
-1.0000
0.0000
1.0000
2.0000

auto-correlation of xn:   size: 5 by 1
0.0000
2.0000
5.0000
2.0000
0.0000

biased auto-correlation of xn:   size: 5 by 1
0.0000
0.6667
1.6667
0.6667
0.0000

unbiased auto-correlation of xn:   size: 5 by 1
0.0000
1.0000
1.6667
1.0000
0.0000

cross-correlation of xn and yn:   size: 9 by 1
0.0000
2.0000
5.0000
2.0000
-1.0000
-4.0000
-4.0000
0.0000
0.0000

cross-correlation of yn and xn:   size: 9 by 1
0.0000
0.0000
-4.0000
-4.0000
-1.0000
2.0000
5.0000
2.0000
0.0000

biased cross-correlation of xn and yn:   size: 9 by 1
0.0000
0.4000
1.0000
0.4000
-0.2000
-0.8000
-0.8000
0.0000
0.0000

biased cross-correlation of yn and xn:   size: 9 by 1
0.0000
0.0000
-0.8000
-0.8000
-0.2000
0.4000
1.0000
0.4000
0.0000

unbiased cross-correlation of xn and yn:   size: 9 by 1
0.0000
1.0000
1.6667
0.5000
-0.2000
-1.0000
-1.3333
0.0000
0.0000

unbiased cross-correlation of yn and xn:   size: 9 by 1
0.0000
0.0000
-1.3333
-1.0000
-0.2000
0.5000
1.6667
1.0000
0.0000

fast auto-correlation of xn:   size: 5 by 1
0.0000
2.0000
5.0000
2.0000
-0.0000

fast biased auto-correlation of xn:   size: 5 by 1
0.0000
0.6667
1.6667
0.6667
-0.0000

fast unbiased auto-correlation of xn:   size: 5 by 1
0.0000
1.0000
1.6667
1.0000
-0.0000

fast cross-correlation of xn and yn:   size: 9 by 1
-0.0000
2.0000
5.0000
2.0000
-1.0000
-4.0000
-4.0000
-0.0000
-0.0000

fast cross-correlation of yn and xn:   size: 9 by 1
-0.0000
-0.0000
-4.0000
-4.0000
-1.0000
2.0000
5.0000
2.0000
0.0000

fast biased cross-correlation of xn and yn:   size: 9 by 1
-0.0000
0.4000
1.0000
0.4000
-0.2000
-0.8000
-0.8000
-0.0000
-0.0000

fast biased cross-correlation of yn and xn:   size: 9 by 1
-0.0000
-0.0000
-0.8000
-0.8000
-0.2000
0.4000
1.0000
0.4000
0.0000

fast unbiased cross-correlation of xn and yn:   size: 9 by 1
-0.0000
1.0000
1.6667
0.5000
-0.2000
-1.0000
-1.3333
-0.0000
-0.0000

fast unbiased cross-correlation of yn and xn:   size: 9 by 1
-0.0000
-0.0000
-1.3333
-1.0000
-0.2000
0.5000
1.6667
1.0000
0.0000


Process returned 0 (0x0)   execution time : 0.187 s
Press any key to continue.

转载于:https://my.oschina.net/zmjerry/blog/4082