This forum has been archived. All content is frozen. Please use KDE Discuss instead.

Cross Correlation function with EigenFFT

Tags: None
(comma "," separated)
kerdemdemir
Registered Member
Posts
2
Karma
0
Hi Everybody ,

I have needed cross correlation function which I wanted to use to correlate 2 sounds. It took a little time searching it so I implemented it my self with EigenFFT.

It takes two vectors and return the index which they correlate most and their correlation value.

I tested it with two gaussian pulses which I added in the code and it worked in my tests .

I wanted to share it , I know it is something easy but took my time . I hope it will be usefull to somebody.

Best Regards

Code: Select all
#include <QCoreApplication>
#include <Eigen/Dense>
#include <unsupported/Eigen/FFT>
#include <deque>
#include <vector>
#include <iostream>

#define PI 3.147

static const double samplesPerSec = 44000;
static double delayTest = -100;


std::vector<double> xCorrTest1;
std::vector<double> xCorrTest2;

void createTestGaussianVector()
{
    //Will create 2 Gaussian Pulses, one of the will have delay according the value of delayTest
    if (delayTest > 0)
    {
        for (int i = 0; i < delayTest; i++)
        {
            xCorrTest2.push_back(0);
        }
    }
    else
    {
        delayTest = delayTest * (-1);
        for (int i = 0; i < delayTest; i++)
        {
            xCorrTest1.push_back(0);
        }
    }
    for (double i = 0; i < (samplesPerSec)*0.15; i++)
    {
        double ts = 1.0 / (double)samplesPerSec;
        double to = 1000;
        double t = (double)((i - ((double)samplesPerSec * 0.075)) * 2.0 * PI * ts) / to;
        xCorrTest1.push_back( 10 * cos(2.0*PI*samplesPerSec*t)
                               * exp(-1.0*(pow(t,2)/pow((ts),2))));
        xCorrTest2.push_back( 10 * cos(2.0*PI*samplesPerSec*t)
                               * exp(-1.0*(pow(t,2)/pow((ts),2))));
    }
}

std::pair<double, double> eigenCrossCorrelation(std::vector<double>& xCorrInputVecFirst
                                                , std::vector<double>& xCorrInputVecSecond)
{
    Eigen::FFT<double> fft;
    int N = std::max(xCorrInputVecFirst.size(), xCorrInputVecSecond.size());

    //Compute the FFT size as the "next power of 2" of the input vector's length (max)
    int b = ceil(log2(2.0 * N - 1));
    int fftsize = pow(2,b);
    int end = fftsize - 1;
    int maxlag = N - 1;
    size_t firstSize = xCorrInputVecFirst.size();
    size_t secondSize = xCorrInputVecSecond.size();

    //Zero Padd
    for (int i = xCorrInputVecFirst.size(); i < fftsize; i++)
    {
        xCorrInputVecFirst.push_back(0);
    }

    for (int i = xCorrInputVecSecond.size(); i < fftsize; i++)
    {
        xCorrInputVecSecond.push_back(0);
    }

    std::vector<std::complex<double> > freqvec;
    std::vector<std::complex<double> > freqvec2;

    //FFT for freq domain to both vectors
    fft.fwd( freqvec,xCorrInputVecFirst);
    fft.fwd( freqvec2,xCorrInputVecSecond);

    //Main step of cross corr
    for (int i = 0; i < fftsize; i++)
    {
        freqvec[i] = freqvec[i] * std::conj(freqvec2[i]);
    }

    std::vector<double> result;
    fft.inv(result, freqvec);

    //Will get rid of extra zero padding and move minus lags to beginning without copy
    std::vector<double> result2(std::make_move_iterator(result.begin() + end - maxlag + 1),
                                std::make_move_iterator(result.end()));

    result2.insert(result2.end(), make_move_iterator(result.begin())
                   , make_move_iterator(result.begin()+maxlag));


    auto minMaxRange = std::minmax_element(result2.begin(),
                        result2.end());

    //Will take back the changes which made in input vector
    if (xCorrInputVecFirst.size() != firstSize)
        xCorrInputVecFirst.resize(firstSize);
    if (xCorrInputVecSecond.size() != secondSize)
        xCorrInputVecSecond.resize(secondSize);


    //Return val
    auto resultIndex = ((minMaxRange.second - result2.begin()) - N + 1);
    std::cout << "Max element at: "  << resultIndex << std::endl;
    auto maxValue = result[minMaxRange.second - result.begin()];
    return std::make_pair (resultIndex, maxValue);
}

int main(int argc, char *argv[])
{
    QCoreApplication a(argc, argv);

    createTestGaussianVector();
    eigenCrossCorrelation(xCorrTest1, xCorrTest2);


    return a.exec();
}

Last edited by kerdemdemir on Sat Nov 30, 2013 10:49 pm, edited 1 time in total.
User avatar
ggael
Moderator
Posts
3447
Karma
19
OS
Hi,

thanks for sharing!


Bookmarks



Who is online

Registered users: Bing [Bot], Google [Bot], Yahoo [Bot]