Added a FFT based correlation class

sdrbase/CMakeLists.txt

@@ -28,6 +28,7 @@ set(sdrbase_SOURCES
     dsp/dspengine.cpp
     dsp/dspdevicesourceengine.cpp
     dsp/dspdevicesinkengine.cpp
+    dsp/fftcorr.cpp
     dsp/fftengine.cpp
     dsp/fftfilt.cxx
     dsp/fftwindow.cpp
@@ -118,6 +119,7 @@ set(sdrbase_HEADERS
     dsp/dspdevicesourceengine.h
     dsp/dspdevicesinkengine.h
     dsp/dsptypes.h
+    dsp/fftcorr.h
     dsp/fftengine.h
     dsp/fftfilt.h
     dsp/fftwengine.h

sdrbase/dsp/fftcorr.cpp

@@ -0,0 +1,112 @@
+///////////////////////////////////////////////////////////////////////////////////
+// Copyright (C) 2018 F4EXB                                                      //
+// written by Edouard Griffiths                                                  //
+//                                                                               //
+// FFT based cross correlation                                                   //
+//                                                                               //
+// See: http://liquidsdr.org/blog/pll-howto/                                     //
+// Fixed filter registers saturation                                             //
+// Added order for PSK locking. This brilliant idea actually comes from this     //
+// post: https://www.dsprelated.com/showthread/comp.dsp/36356-1.php              //
+//                                                                               //
+// 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 as version 3 of the License, or                  //
+//                                                                               //
+// 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 V3 for more details.                               //
+//                                                                               //
+// You should have received a copy of the GNU General Public License             //
+// along with this program. If not, see <http://www.gnu.org/licenses/>.          //
+///////////////////////////////////////////////////////////////////////////////////
+
+#include <algorithm>
+#include "fftcorr.h"
+
+void fftcorr::init_fft()
+{
+    flen2   = flen >> 1;
+    fftA = new g_fft<float>(flen);
+    fftB = new g_fft<float>(flen);
+
+    dataA   = new cmplx[flen];
+    dataB   = new cmplx[flen];
+    dataBj  = new cmplx[flen];
+    dataP   = new cmplx[flen];
+    output  = new cmplx[flen2];
+    ovlbuf  = new cmplx[flen2];
+
+    std::fill(dataA, dataA+flen, 0);
+    std::fill(dataB, dataB+flen, 0);
+    std::fill(output, output+flen2, 0);
+    std::fill(ovlbuf, ovlbuf+flen2, 0);
+
+    inptrA = 0;
+    inptrB = 0;
+}
+
+fftcorr::fftcorr(int len) : flen(len), flen2(len>>1)
+{
+    init_fft();
+}
+
+fftcorr::~fftcorr()
+{
+    delete fftA;
+    delete fftB;
+    delete[] dataA;
+    delete[] dataB;
+    delete[] dataBj;
+    delete[] dataP;
+    delete[] output;
+    delete[] ovlbuf;
+}
+
+int fftcorr::run(const cmplx& inA, const cmplx* inB, cmplx **out)
+{
+    dataA[inptrA++] = inA;
+
+    if (inB) {
+        dataB[inptrB++] = *inB;
+    }
+
+    if (inptrA < flen2) {
+        return 0;
+    }
+
+    fftA->ComplexFFT(dataA);
+
+    if (inB) {
+        fftB->ComplexFFT(dataB);
+    }
+
+    if (inB) {
+        std::transform(dataB, dataB+flen, dataBj, [](const cmplx& c) -> cmplx { return std::conj(c); });
+    } else {
+        std::transform(dataA, dataA+flen, dataBj, [](const cmplx& c) -> cmplx { return std::conj(c); });
+    }
+
+    std::transform(dataA, dataA+flen, dataBj, dataP, [](const cmplx& a, const cmplx& b) -> cmplx { return a*b; });
+
+    fftA->InverseComplexFFT(dataP);
+
+    for (int i = 0; i < flen2; i++)
+    {
+        output[i] = ovlbuf[i] + dataP[i];
+        ovlbuf[i] = dataP[flen2 + i];
+    }
+
+    std::fill(dataA, dataA+flen, 0);
+    inptrA = 0;
+
+    if (inB)
+    {
+        std::fill(dataB, dataB+flen, 0);
+        inptrB = 0;
+    }
+
+    *out = output;
+    return flen2;
+}

sdrbase/dsp/fftcorr.h

@@ -0,0 +1,57 @@
+///////////////////////////////////////////////////////////////////////////////////
+// Copyright (C) 2018 F4EXB                                                      //
+// written by Edouard Griffiths                                                  //
+//                                                                               //
+// FFT based cross correlation                                                   //
+//                                                                               //
+// See: http://liquidsdr.org/blog/pll-howto/                                     //
+// Fixed filter registers saturation                                             //
+// Added order for PSK locking. This brilliant idea actually comes from this     //
+// post: https://www.dsprelated.com/showthread/comp.dsp/36356-1.php              //
+//                                                                               //
+// 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 as version 3 of the License, or                  //
+//                                                                               //
+// 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 V3 for more details.                               //
+//                                                                               //
+// You should have received a copy of the GNU General Public License             //
+// along with this program. If not, see <http://www.gnu.org/licenses/>.          //
+///////////////////////////////////////////////////////////////////////////////////
+
+#ifndef SDRBASE_DSP_FFTCORR_H_
+#define SDRBASE_DSP_FFTCORR_H_
+
+#include <complex>
+#include "gfft.h"
+#include "export.h"
+
+class SDRBASE_API fftcorr {
+public:
+    typedef std::complex<float> cmplx;
+    fftcorr(int len);
+    ~fftcorr();
+
+    int run(const cmplx& inA, const cmplx* inB, cmplx **out); //!< if inB = 0 then run auto-correlation
+
+private:
+    void init_fft();
+    int flen;  //!< FFT length
+    int flen2; //!< half FFT length
+    g_fft<float> *fftA;
+    g_fft<float> *fftB;
+    cmplx *dataA;  // from A input
+    cmplx *dataB;  // from B input
+    cmplx *dataBj; // conjugate of B
+    cmplx *dataP;  // product of A with conjugate of B
+    cmplx *ovlbuf;
+    cmplx *output;
+    int inptrA;
+    int inptrB;
+};
+
+
+#endif /* SDRBASE_DSP_FFTCORR_H_ */