From d673278f86c21c4579346c48a94b7a6d51f5b21e Mon Sep 17 00:00:00 2001 From: f4exb Date: Sun, 20 May 2018 10:42:14 +0200 Subject: [PATCH] Added a FFT based correlation class --- sdrbase/CMakeLists.txt | 2 + sdrbase/dsp/fftcorr.cpp | 112 ++++++++++++++++++++++++++++++++++++++++ sdrbase/dsp/fftcorr.h | 57 ++++++++++++++++++++ 3 files changed, 171 insertions(+) create mode 100644 sdrbase/dsp/fftcorr.cpp create mode 100644 sdrbase/dsp/fftcorr.h diff --git a/sdrbase/CMakeLists.txt b/sdrbase/CMakeLists.txt index 5bacdd135..da49b3469 100644 --- a/sdrbase/CMakeLists.txt +++ b/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 diff --git a/sdrbase/dsp/fftcorr.cpp b/sdrbase/dsp/fftcorr.cpp new file mode 100644 index 000000000..22fa00ba4 --- /dev/null +++ b/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 . // +/////////////////////////////////////////////////////////////////////////////////// + +#include +#include "fftcorr.h" + +void fftcorr::init_fft() +{ + flen2 = flen >> 1; + fftA = new g_fft(flen); + fftB = new g_fft(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; +} diff --git a/sdrbase/dsp/fftcorr.h b/sdrbase/dsp/fftcorr.h new file mode 100644 index 000000000..cd7777340 --- /dev/null +++ b/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 . // +/////////////////////////////////////////////////////////////////////////////////// + +#ifndef SDRBASE_DSP_FFTCORR_H_ +#define SDRBASE_DSP_FFTCORR_H_ + +#include +#include "gfft.h" +#include "export.h" + +class SDRBASE_API fftcorr { +public: + typedef std::complex 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 *fftA; + g_fft *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_ */