Channel Analyzer: use FFTW (or Kiss) engine with Hamming window for auto correlation

This commit is contained in:
f4exb 2019-08-18 01:58:09 +02:00
parent 7b8319c79f
commit e12f4ed71b
5 changed files with 186 additions and 4 deletions

View File

@ -52,7 +52,7 @@ ChannelAnalyzer::ChannelAnalyzer(DeviceAPI *deviceAPI) :
SSBFilter = new fftfilt(m_settings.m_lowCutoff / m_inputSampleRate, m_settings.m_bandwidth / m_inputSampleRate, ssbFftLen);
DSBFilter = new fftfilt(m_settings.m_bandwidth / m_inputSampleRate, 2*ssbFftLen);
RRCFilter = new fftfilt(m_settings.m_bandwidth / m_inputSampleRate, 2*ssbFftLen);
m_corr = new fftcorr(8*ssbFftLen); // 8k for 4k effective samples
m_corr = new fftcorr2(8*ssbFftLen); // 8k for 4k effective samples
m_pll.computeCoefficients(0.002f, 0.5f, 10.0f); // bandwidth, damping factor, loop gain
applyChannelSettings(m_inputSampleRate, m_inputFrequencyOffset, true);
@ -73,6 +73,7 @@ ChannelAnalyzer::~ChannelAnalyzer()
delete SSBFilter;
delete DSBFilter;
delete RRCFilter;
delete m_corr;
}
void ChannelAnalyzer::feed(const SampleVector::const_iterator& begin, const SampleVector::const_iterator& end, bool positiveOnly)

View File

@ -25,7 +25,7 @@
#include "channel/channelapi.h"
#include "dsp/interpolator.h"
#include "dsp/ncof.h"
#include "dsp/fftcorr.h"
#include "dsp/fftcorr2.h"
#include "dsp/fftfilt.h"
#include "dsp/phaselockcomplex.h"
#include "dsp/freqlockcomplex.h"
@ -181,7 +181,7 @@ private:
fftfilt* SSBFilter;
fftfilt* DSBFilter;
fftfilt* RRCFilter;
fftcorr* m_corr;
fftcorr2* m_corr;
BasebandSampleSink* m_sampleSink;
SampleVector m_sampleBuffer;
@ -209,7 +209,8 @@ private:
break;
case ChannelAnalyzerSettings::InputAutoCorr:
{
std::complex<float> a = m_corr->run(s/(SDR_RX_SCALEF/768.0f), 0);
//std::complex<float> a = m_corr->run(s/(SDR_RX_SCALEF/768.0f), 0);
std::complex<float> a = m_corr->run(s/SDR_RX_SCALEF, 0);
if (m_settings.m_ssb & !m_usb) { // invert spectrum for LSB
m_sampleBuffer.push_back(Sample(a.imag(), a.real()));

View File

@ -77,6 +77,7 @@ set(sdrbase_SOURCES
dsp/dspdevicesinkengine.cpp
dsp/dspdevicemimoengine.cpp
dsp/fftcorr.cpp
dsp/fftcorr2.cpp
dsp/fftengine.cpp
dsp/fftfilt.cpp
dsp/fftwindow.cpp
@ -193,6 +194,7 @@ set(sdrbase_HEADERS
dsp/dspdevicemimoengine.h
dsp/dsptypes.h
dsp/fftcorr.h
dsp/fftcorr2.h
dsp/fftengine.h
dsp/fftfilt.h
dsp/fftwengine.h

120
sdrbase/dsp/fftcorr2.cpp Normal file
View File

@ -0,0 +1,120 @@
///////////////////////////////////////////////////////////////////////////////////
// Copyright (C) 2019 F4EXB //
// written by Edouard Griffiths //
// //
// FFT based cross correlation. Uses FFTW/Kiss engine. //
// //
// 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 //
// (at your option) any later version. //
// //
// 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 "dsp/fftengine.h"
#include "fftcorr2.h"
void fftcorr2::init_fft()
{
fftA->configure(flen, false);
fftB->configure(flen, false);
fftInvA->configure(flen, true);
m_window.create(FFTWindow::Hamming, flen);
dataA = new cmplx[flen];
dataB = new cmplx[flen];
dataBj = new cmplx[flen];
dataP = new cmplx[flen];
std::fill(dataA, dataA+flen, 0);
std::fill(dataB, dataB+flen, 0);
inptrA = 0;
inptrB = 0;
outptr = 0;
}
fftcorr2::fftcorr2(int len) :
flen(len),
flen2(len>>1),
fftA(FFTEngine::create()),
fftB(FFTEngine::create()),
fftInvA(FFTEngine::create())
{
init_fft();
}
fftcorr2::~fftcorr2()
{
delete fftA;
delete fftB;
delete[] dataA;
delete[] dataB;
delete[] dataBj;
delete[] dataP;
}
int fftcorr2::run(const cmplx& inA, const cmplx* inB, cmplx **out)
{
dataA[inptrA++] = inA;
if (inB) {
dataB[inptrB++] = *inB;
}
if (inptrA < flen2) {
return 0;
}
m_window.apply(dataA, fftA->in());
fftA->transform();
if (inB)
{
m_window.apply(dataB, fftB->in());
fftB->transform();
}
if (inB) {
std::transform(fftB->out(), fftB->out()+flen, dataBj, [](const cmplx& c) -> cmplx { return std::conj(c); });
} else {
std::transform(fftA->out(), fftA->out()+flen, dataBj, [](const cmplx& c) -> cmplx { return std::conj(c); });
}
std::transform(fftA->out(), fftA->out()+flen, dataBj, fftInvA->in(), [](const cmplx& a, const cmplx& b) -> cmplx { return a*b; });
fftInvA->transform();
std::copy(fftInvA->out(), fftInvA->out()+flen, dataP);
std::fill(dataA, dataA+flen, 0);
inptrA = 0;
if (inB)
{
std::fill(dataB, dataB+flen, 0);
inptrB = 0;
}
*out = dataP;
return flen2;
}
const fftcorr2::cmplx& fftcorr2::run(const cmplx& inA, const cmplx* inB)
{
cmplx *dummy;
if (run(inA, inB, &dummy)) {
outptr = 0;
}
return dataP[outptr++];
}

58
sdrbase/dsp/fftcorr2.h Normal file
View File

@ -0,0 +1,58 @@
///////////////////////////////////////////////////////////////////////////////////
// Copyright (C) 2019 F4EXB //
// written by Edouard Griffiths //
// //
// FFT based cross correlation. Uses FFTW/Kiss engine. //
// //
// 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 //
// (at your option) any later version. //
// //
// 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_FFTCORR2_H_
#define SDRBASE_DSP_FFTCORR2_H_
#include <complex>
#include "dsp/fftwindow.h"
#include "export.h"
class FFTEngine;
class SDRBASE_API fftcorr2 {
public:
typedef std::complex<float> cmplx;
fftcorr2(int len);
~fftcorr2();
int run(const cmplx& inA, const cmplx* inB, cmplx **out); //!< if inB = 0 then run auto-correlation
const cmplx& run(const cmplx& inA, const cmplx* inB);
private:
void init_fft();
int flen; //!< FFT length
int flen2; //!< half FFT length
FFTEngine *fftA;
FFTEngine *fftB;
FFTEngine *fftInvA;
FFTWindow m_window;
cmplx *dataA; // from A input
cmplx *dataB; // from B input
cmplx *dataBj; // conjugate of B
cmplx *dataP; // product of A with conjugate of B
int inptrA;
int inptrB;
int outptr;
};
#endif /* SDRBASE_DSP_FFTCORR2_H_ */