xcorr_cpp

对两个等长的复数数组做xcorr，效果与matlab xcorr函数保持一致

存在两个.mat文件，文件中存储着N*1个复数，对两个数组进行xcorr。首先使用matlab的mat处理接口读取两个文件的数据，然后求xcorr，最后绘图或保存数据。

// main.cpp
#include <QApplication>
#include <QDebug>
#include <QFile>
#include <QTextStream>
#include <QVector>
#include <qmath.h>
#include "mat.h"
#include "matrix.h"
#include "qcustomPlot.h"

#define PLOT
// #define WRITE
// #define TEST

/*
if i = 0 , x[0]*y[n-1] 1
if i = 1 , x[0]*y[n-2]+x[1]*y[n-1] 2
if i = 2 , x[0]*y[n-3]+x[1]*y[n-2]+x[2]*y[n-1] 3
if i = 3 , x[0]*y[n-4]+x[1]*y[n-3]+x[2]*y[n-2]+x[3]*y[n-1] 4
...
if i = n-1 , x[0]*y[0]+x[1]*y[1]+x[2]*y[2]+...+x[n-1]*y[n-1] n
if i = n   , x[1]*y[0]+x[2]*y[1]+x[3]*y[2]+...+x[n-1]*y[n-2] n-1
if i = n+1 , x[2]*y[0]+x[3]*y[1]+x[4]*y[2]+...+x[n-1]*y[n-3] n-2
...
if i = 2n-2 , x[n-1]*y[0] 1
*/

mxArray *xcorr(const mxArray *x, const mxArray *y)
{
    auto *p1 = mxGetComplexDoubles(x);
    auto *p2 = mxGetComplexDoubles(y);
    size_t m = mxGetM(x);
    mxArray *out = mxCreateDoubleMatrix(2 * m - 1, 1, mxCOMPLEX);
    auto *p3 = mxGetComplexDoubles(out);
    for (size_t i = 0; i < 2 * m - 1; ++i)
    {
        p3[i].real = 0;
        p3[i].imag = 0;
        size_t start = i <= m - 1 ? 0 : i - m + 1;
        size_t end = i <= m - 1 ? i : m - 1;
        for (size_t j = start; j <= end; ++j)
        {
            int k = m - 1 - i + j;
            p3[i].real += p1[j].real * p2[k].real + p1[j].imag * p2[k].imag;
            p3[i].imag += p1[j].imag * p2[k].real - p1[j].real * p2[k].imag;
        }
    }
    return out;
}
mxArray *readData(const char *file, const char *varname)
{
    MATFile *pmat = matOpen(file, "r");
    mxArray *pa = nullptr;
    if (pmat)
    {
        pa = matGetVariable(pmat, varname);
        matClose(pmat);
    }
    return pa;
}

int main(int argc, char *argv[])
{
    QApplication app(argc, argv);

    auto chanel_1 = readData("path","name"); // 根据自己的路径修改path、name
    auto chanel_2 = readData("path","name");
#ifdef PLOT
    auto out = xcorr(chanel_1, chanel_2);
    size_t m = mxGetM(out);
    QVector<double> plot_data(m), x_data(m);
    for (size_t i = 0; i < m; ++i)
    {
        x_data[i] = i;
        plot_data[i] = qSqrt(mxGetComplexDoubles(out)[i].real * mxGetComplexDoubles(out)[i].real + mxGetComplexDoubles(out)[i].imag * mxGetComplexDoubles(out)[i].imag);
    }

    QCustomPlot plot;
    // plot.setOpenGl(true);
    plot.addGraph();
    plot.graph(0)->setData(x_data, plot_data);
    plot.graph(0)->setPen(QPen(Qt::red));
    plot.rescaleAxes();
    plot.replot();
    plot.show();
    plot.setInteractions(QCP::iRangeDrag | QCP::iRangeZoom | QCP::iSelectPlottables);
#endif // PLOT

#ifdef WRITE
    auto *out = xcorr(chanel_1, chanel_2);
    size_t m = mxGetM(out);
    auto *p = mxGetComplexDoubles(out);
    QFile file("xcorr_data.txt");
    if (file.open(QIODevice::WriteOnly | QIODevice::Text))
    {
        QTextStream out(&file);
        for (size_t i = 0; i < m; ++i)
        {
            out << p[i].real << " " << p[i].imag << "\n";
        }
    }
    qDebug() << "write done";
    file.close();
    mxFree(p);
#endif // WRITE

#ifdef TEST
    auto *p1 = mxGetComplexDoubles(chanel_1);
    auto *p2 = mxGetComplexDoubles(chanel_2);
    size_t m = mxGetM(chanel_1);
    QVector<double> data_real(m), data_imag(m);
    for (size_t i = 0; i < m; ++i)
    {
        data_real[i] = p1[i].real * p2[i].real + p1[i].imag * p2[i].imag;
        data_imag[i] = p1[i].imag * p2[i].real - p1[i].real * p2[i].imag;
    }
#endif // TEST

    return app.exec();
}