2018-05-20 08:42:14 +00:00
|
|
|
///////////////////////////////////////////////////////////////////////////////////
|
2019-10-27 21:25:05 +00:00
|
|
|
// Copyright (C) 2019 F4EXB //
|
2018-05-20 08:42:14 +00:00
|
|
|
// written by Edouard Griffiths //
|
|
|
|
// //
|
2019-10-27 21:25:05 +00:00
|
|
|
// FFT based cross correlation. Uses FFTW/Kiss engine. //
|
2018-05-20 08:42:14 +00:00
|
|
|
// //
|
|
|
|
// 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 //
|
2019-04-11 12:32:15 +00:00
|
|
|
// (at your option) any later version. //
|
2018-05-20 08:42:14 +00:00
|
|
|
// //
|
|
|
|
// 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>
|
2019-10-27 21:25:05 +00:00
|
|
|
|
|
|
|
#include "dsp/fftengine.h"
|
2018-05-20 08:42:14 +00:00
|
|
|
#include "fftcorr.h"
|
|
|
|
|
|
|
|
void fftcorr::init_fft()
|
|
|
|
{
|
2019-10-27 21:25:05 +00:00
|
|
|
fftA->configure(flen, false);
|
|
|
|
fftB->configure(flen, false);
|
|
|
|
fftInvA->configure(flen, true);
|
|
|
|
m_window.create(FFTWindow::Hanning, flen);
|
2018-05-20 08:42:14 +00:00
|
|
|
|
2019-10-27 21:25:05 +00:00
|
|
|
dataA = new cmplx[flen];
|
|
|
|
dataB = new cmplx[flen];
|
|
|
|
dataBj = new cmplx[flen];
|
|
|
|
dataP = new cmplx[flen];
|
2018-05-20 08:42:14 +00:00
|
|
|
|
|
|
|
std::fill(dataA, dataA+flen, 0);
|
|
|
|
std::fill(dataB, dataB+flen, 0);
|
|
|
|
|
|
|
|
inptrA = 0;
|
|
|
|
inptrB = 0;
|
2018-05-20 16:17:53 +00:00
|
|
|
outptr = 0;
|
2018-05-20 08:42:14 +00:00
|
|
|
}
|
|
|
|
|
2019-10-27 21:25:05 +00:00
|
|
|
fftcorr::fftcorr(int len) :
|
|
|
|
flen(len),
|
|
|
|
flen2(len>>1),
|
|
|
|
fftA(FFTEngine::create()),
|
|
|
|
fftB(FFTEngine::create()),
|
|
|
|
fftInvA(FFTEngine::create())
|
2018-05-20 08:42:14 +00:00
|
|
|
{
|
|
|
|
init_fft();
|
|
|
|
}
|
|
|
|
|
|
|
|
fftcorr::~fftcorr()
|
|
|
|
{
|
|
|
|
delete fftA;
|
|
|
|
delete fftB;
|
|
|
|
delete[] dataA;
|
|
|
|
delete[] dataB;
|
|
|
|
delete[] dataBj;
|
|
|
|
delete[] dataP;
|
|
|
|
}
|
|
|
|
|
|
|
|
int fftcorr::run(const cmplx& inA, const cmplx* inB, cmplx **out)
|
|
|
|
{
|
|
|
|
dataA[inptrA++] = inA;
|
|
|
|
|
|
|
|
if (inB) {
|
|
|
|
dataB[inptrB++] = *inB;
|
|
|
|
}
|
|
|
|
|
2018-05-20 18:23:41 +00:00
|
|
|
if (inptrA < flen2) {
|
2018-05-20 08:42:14 +00:00
|
|
|
return 0;
|
|
|
|
}
|
|
|
|
|
2019-10-27 21:25:05 +00:00
|
|
|
m_window.apply(dataA, fftA->in());
|
|
|
|
fftA->transform();
|
2018-05-20 08:42:14 +00:00
|
|
|
|
2019-10-27 21:25:05 +00:00
|
|
|
if (inB)
|
|
|
|
{
|
|
|
|
m_window.apply(dataB, fftB->in());
|
|
|
|
fftB->transform();
|
2018-05-20 08:42:14 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
if (inB) {
|
2019-10-27 21:25:05 +00:00
|
|
|
std::transform(fftB->out(), fftB->out()+flen, dataBj, [](const cmplx& c) -> cmplx { return std::conj(c); });
|
2018-05-20 08:42:14 +00:00
|
|
|
} else {
|
2019-10-27 21:25:05 +00:00
|
|
|
std::transform(fftA->out(), fftA->out()+flen, dataBj, [](const cmplx& c) -> cmplx { return std::conj(c); });
|
2018-05-20 08:42:14 +00:00
|
|
|
}
|
|
|
|
|
2019-10-27 21:25:05 +00:00
|
|
|
std::transform(fftA->out(), fftA->out()+flen, dataBj, fftInvA->in(), [](const cmplx& a, const cmplx& b) -> cmplx { return a*b; });
|
2018-05-20 18:23:41 +00:00
|
|
|
|
2019-10-27 21:25:05 +00:00
|
|
|
fftInvA->transform();
|
|
|
|
std::copy(fftInvA->out(), fftInvA->out()+flen, dataP);
|
2018-05-20 08:42:14 +00:00
|
|
|
|
|
|
|
std::fill(dataA, dataA+flen, 0);
|
|
|
|
inptrA = 0;
|
|
|
|
|
|
|
|
if (inB)
|
|
|
|
{
|
|
|
|
std::fill(dataB, dataB+flen, 0);
|
|
|
|
inptrB = 0;
|
|
|
|
}
|
|
|
|
|
2018-05-20 17:41:36 +00:00
|
|
|
*out = dataP;
|
2018-05-20 18:23:41 +00:00
|
|
|
return flen2;
|
2018-05-20 08:42:14 +00:00
|
|
|
}
|
2018-05-20 16:17:53 +00:00
|
|
|
|
|
|
|
const fftcorr::cmplx& fftcorr::run(const cmplx& inA, const cmplx* inB)
|
|
|
|
{
|
|
|
|
cmplx *dummy;
|
|
|
|
|
|
|
|
if (run(inA, inB, &dummy)) {
|
|
|
|
outptr = 0;
|
|
|
|
}
|
|
|
|
|
2018-05-20 17:41:36 +00:00
|
|
|
return dataP[outptr++];
|
2018-05-20 16:17:53 +00:00
|
|
|
}
|