From a5e4e28f067d1d86a127146ede49ded4ed3b4a0d Mon Sep 17 00:00:00 2001
From: Rob van Nieuwpoort <nieuwpoort@astron.nl>
Date: Wed, 27 Jun 2012 10:00:34 +0000
Subject: [PATCH] Task #1611: working on online flagger that flags in frequency
even if there are no channels
---
.gitattributes | 2 +
RTCP/CNProc/src/CMakeLists.txt | 1 +
RTCP/CNProc/src/CN_Processing.cc | 26 ++
RTCP/CNProc/src/CN_Processing.h | 3 +
RTCP/CNProc/src/InversePPF.cc | 17 +-
RTCP/CNProc/src/InversePPF.h | 6 +-
RTCP/CNProc/src/PPF.cc | 2 +-
.../src/PreCorrelationNoChannelsFlagger.cc | 270 ++++++++++++++++++
.../src/PreCorrelationNoChannelsFlagger.h | 63 ++++
RTCP/CNProc/test/FlaggerTest.parset | 14 +-
RTCP/CNProc/test/RFI-test.parset | 4 +-
RTCP/CNProc/test/tInversePPF.cc | 9 +-
RTCP/Interface/include/Interface/Parset.h | 6 +
13 files changed, 397 insertions(+), 26 deletions(-)
create mode 100644 RTCP/CNProc/src/PreCorrelationNoChannelsFlagger.cc
create mode 100644 RTCP/CNProc/src/PreCorrelationNoChannelsFlagger.h
diff --git a/.gitattributes b/.gitattributes
index 86e1ca97270..5251397a6a0 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -3421,6 +3421,8 @@ RTCP/CNProc/src/PostCorrelationFlagger.cc -text
RTCP/CNProc/src/PostCorrelationFlagger.h -text
RTCP/CNProc/src/PreCorrelationFlagger.cc -text
RTCP/CNProc/src/PreCorrelationFlagger.h -text
+RTCP/CNProc/src/PreCorrelationNoChannelsFlagger.cc -text
+RTCP/CNProc/src/PreCorrelationNoChannelsFlagger.h -text
RTCP/CNProc/src/Ring.h -text
RTCP/CNProc/src/Stokes.cc -text
RTCP/CNProc/src/Stokes.h -text
diff --git a/RTCP/CNProc/src/CMakeLists.txt b/RTCP/CNProc/src/CMakeLists.txt
index 1dcdbfdc2da..3a868609839 100644
--- a/RTCP/CNProc/src/CMakeLists.txt
+++ b/RTCP/CNProc/src/CMakeLists.txt
@@ -29,6 +29,7 @@ set(cnproc_LIB_SRCS
PPF.cc
Flagger.cc
PreCorrelationFlagger.cc
+ PreCorrelationNoChannelsFlagger.cc
PostCorrelationFlagger.cc
Stokes.cc)
diff --git a/RTCP/CNProc/src/CN_Processing.cc b/RTCP/CNProc/src/CN_Processing.cc
index 82d68a74877..b8cac2aa629 100644
--- a/RTCP/CNProc/src/CN_Processing.cc
+++ b/RTCP/CNProc/src/CN_Processing.cc
@@ -208,6 +208,12 @@ template <typename SAMPLE_TYPE> CN_Processing<SAMPLE_TYPE>::CN_Processing(const
LOG_DEBUG_STR("Online PreCorrelation flagger enabled");
}
+ if (parset.onlineFlagging() && parset.onlinePreCorrelationNoChannelsFlagging()) {
+ itsPreCorrelationNoChannelsFlagger = new PreCorrelationNoChannelsFlagger(parset, itsNrStations, itsNrSubbands, itsNrChannels, itsNrSamplesPerIntegration);
+ if (LOG_CONDITION)
+ LOG_DEBUG_STR("Online PreCorrelation no channels flagger enabled");
+ }
+
if (parset.outputCorrelatedData()) {
itsCorrelator = new Correlator(itsBeamFormer->getStationMapping(), itsNrChannels, itsNrSamplesPerIntegration);
itsCorrelatedData = (CorrelatedData*)newStreamableData(parset, CORRELATED_DATA, -1, itsBigAllocator);
@@ -726,6 +732,23 @@ template <typename SAMPLE_TYPE> void CN_Processing<SAMPLE_TYPE>::preCorrelationF
}
+template <typename SAMPLE_TYPE> void CN_Processing<SAMPLE_TYPE>::preCorrelationNoChannelsFlagging()
+{
+#if defined HAVE_MPI
+ if (LOG_CONDITION)
+ LOG_DEBUG_STR(itsLogPrefix << "Start pre correlation no channels flagger at t = " << blockAge());
+#endif // HAVE_MPI
+
+ static NSTimer timer("pre correlation no channels flagger", true, true);
+
+ timer.start();
+ computeTimer.start();
+ itsPreCorrelationNoChannelsFlagger->flag(itsFilteredData, *itsCurrentSubband);
+ computeTimer.stop();
+ timer.stop();
+}
+
+
template <typename SAMPLE_TYPE> void CN_Processing<SAMPLE_TYPE>::mergeStations()
{
#if defined HAVE_MPI
@@ -930,6 +953,9 @@ template <typename SAMPLE_TYPE> void CN_Processing<SAMPLE_TYPE>::process(unsigne
if (itsPPF != 0)
filter();
+ if (itsPreCorrelationNoChannelsFlagger != 0)
+ preCorrelationNoChannelsFlagging();
+
if (itsPreCorrelationFlagger != 0)
preCorrelationFlagging();
diff --git a/RTCP/CNProc/src/CN_Processing.h b/RTCP/CNProc/src/CN_Processing.h
index 3ab1214816b..a3f10c2a141 100644
--- a/RTCP/CNProc/src/CN_Processing.h
+++ b/RTCP/CNProc/src/CN_Processing.h
@@ -47,6 +47,7 @@
#include <LocationInfo.h>
#include <PPF.h>
#include <PreCorrelationFlagger.h>
+#include <PreCorrelationNoChannelsFlagger.h>
#include <PostCorrelationFlagger.h>
#include <Ring.h>
#include <Stokes.h>
@@ -86,6 +87,7 @@ template <typename SAMPLE_TYPE> class CN_Processing : public CN_Processing_Base
int transposeBeams(unsigned block);
void filter();
void dedisperseAfterBeamForming(unsigned beam, double dm);
+ void preCorrelationNoChannelsFlagging();
void preCorrelationFlagging();
void mergeStations();
void formBeams(unsigned sap, unsigned firstBeam, unsigned nrBeams);
@@ -168,6 +170,7 @@ template <typename SAMPLE_TYPE> class CN_Processing : public CN_Processing_Base
SmartPtr<DedispersionAfterBeamForming> itsDedispersionAfterBeamForming;
SmartPtr<DedispersionBeforeBeamForming> itsDedispersionBeforeBeamForming;
SmartPtr<PreCorrelationFlagger> itsPreCorrelationFlagger;
+ SmartPtr<PreCorrelationNoChannelsFlagger> itsPreCorrelationNoChannelsFlagger;
SmartPtr<PostCorrelationFlagger> itsPostCorrelationFlagger;
SmartPtr<Trigger> itsTrigger;
};
diff --git a/RTCP/CNProc/src/InversePPF.cc b/RTCP/CNProc/src/InversePPF.cc
index 49c5d7c1c53..5bf844a3d33 100644
--- a/RTCP/CNProc/src/InversePPF.cc
+++ b/RTCP/CNProc/src/InversePPF.cc
@@ -76,18 +76,13 @@ void InversePPF::destroyFFT() {
#endif
}
-// Reads itsFftOutData, writes to invertedFilteredData.
-void InversePPF::performFilter(InverseFilteredData& invertedFilteredData, unsigned time, unsigned minorTime) {
- float sample = itsFftOutData[minorTime];
- float result = itsFIRs[minorTime].processNextSample(sample);
- invertedFilteredData.samples[time * itsOnStationFilterSize + minorTime] = result;
-}
-
// Goes from tansposedBeamFormedData to itsFftInData.
void InversePPF::createFFTInput(const TransposedBeamFormedData& transposedBeamFormedData, unsigned time) {
fftInTimer.start();
// First set the unselected subbands to zero.
+ // We have to do this every time, since the input is destroyed by the FFT.
+ // However, the time this takes is very small compared to the time to fill in the real data below.
memset(itsFftInData, 0, itsOnStationFilterSize * sizeof(float));
// Fill input buffer, using "half complex" format.
@@ -98,7 +93,7 @@ void InversePPF::createFFTInput(const TransposedBeamFormedData& transposedBeamFo
fcomplex sample = transposedBeamFormedData.samples[sb][0 /* channel, but there only is 1 now */][time];
itsFftInData[sb] = real(sample);
- itsFftInData[itsOnStationFilterSize - sb] = imag(sample);
+ itsFftInData[itsOnStationFilterSize - sb - 1] = imag(sample);
}
fftInTimer.stop();
@@ -119,11 +114,15 @@ void InversePPF::performInverseFFT() {
fftTimer.stop();
}
+// Reads itsFftOutData, writes to invertedFilteredData.
void InversePPF::performFiltering(InverseFilteredData& invertedFilteredData, unsigned time) {
firTimer.start();
+ unsigned index = time * itsOnStationFilterSize;
for (unsigned minorTime = 0; minorTime < itsOnStationFilterSize; minorTime++) {
- performFilter(invertedFilteredData, time, minorTime);
+ const float sample = itsFftOutData[minorTime];
+ const float result = itsFIRs[minorTime].processNextSample(sample);
+ invertedFilteredData.samples[index++] = result;
}
firTimer.stop();
diff --git a/RTCP/CNProc/src/InversePPF.h b/RTCP/CNProc/src/InversePPF.h
index 00b67065dda..db9ac15e343 100644
--- a/RTCP/CNProc/src/InversePPF.h
+++ b/RTCP/CNProc/src/InversePPF.h
@@ -30,7 +30,7 @@
* Input must be in "half complex" format.
This consists of the non-redundant half of the complex output for a 1d real-input DFT of size n,
- stored as a sequence of n real numbers (double) in the format:
+ stored as a sequence of n real numbers (double) in the format:
r0, r1, r2, ..., rn/2, i(n+1)/2-1, ..., i2, i1
@@ -105,7 +105,6 @@ private:
void initFFT();
void destroyFFT();
- void performFilter(InverseFilteredData& invertedFilteredData, unsigned time, unsigned minorTime);
void createFFTInput(const TransposedBeamFormedData& transposedBeamFormedData, unsigned time);
void performInverseFFT();
void performFiltering(InverseFilteredData& invertedFilteredData, unsigned time);
@@ -127,10 +126,9 @@ private:
unsigned itsNrSubbands;
unsigned itsNrTaps; // 16
unsigned itsNrSamplesPerIntegration;
- unsigned itsOnStationFilterSize; // 1204
+ unsigned itsOnStationFilterSize; // 1024
bool itsVerbose;
-
};
} // namespace RTCP
diff --git a/RTCP/CNProc/src/PPF.cc b/RTCP/CNProc/src/PPF.cc
index 2500173d38d..60054befef8 100644
--- a/RTCP/CNProc/src/PPF.cc
+++ b/RTCP/CNProc/src/PPF.cc
@@ -385,7 +385,7 @@ template <typename SAMPLE_TYPE> void PPF<SAMPLE_TYPE>::bypass(unsigned stat, dou
#if defined PPF_C_IMPLEMENTATION
for (unsigned time = 0; time < itsNrSamplesPerIntegration; time ++) {
for (unsigned pol = 0; pol < NR_POLARIZATIONS; pol ++) {
- if (filteredData->flags[1][stat].test(time)) {
+ if ((itsNrChannels > 1 && filteredData->flags[1][stat].test(time)) || (itsNrChannels == 1 && filteredData->flags[0][stat].test(time))) {
filteredData->samples[0][stat][time][pol] = makefcomplex(0, 0);
} else {
SAMPLE_TYPE currSample = transposedData->samples[stat][time + alignmentShift][pol];
diff --git a/RTCP/CNProc/src/PreCorrelationNoChannelsFlagger.cc b/RTCP/CNProc/src/PreCorrelationNoChannelsFlagger.cc
new file mode 100644
index 00000000000..98124e03c0f
--- /dev/null
+++ b/RTCP/CNProc/src/PreCorrelationNoChannelsFlagger.cc
@@ -0,0 +1,270 @@
+//# Always #include <lofar_config.h> first!
+#include <lofar_config.h>
+
+#include <Common/Timer.h>
+
+#include <PreCorrelationNoChannelsFlagger.h>
+
+// history is kept per subband, as we can get different subbands over time on this compute node.
+// Always flag poth polarizations as a unit.
+
+// FFT followed by an inverse FFT multiplies all samples by N. Thus, we have to divide by N after we are done.
+
+
+/*
+ First, we flag in the time direction, while integrating to imrove signal-to-noise.
+ This was empirically verified to work much better than flagging on the raw data.
+ We can then replace flagged samples with 0s or mean/ median.
+
+ Two options for frequency flagging:
+
+ - integrate until we have FFTSize samples, so we improve signal-to-noise
+ - do FFT
+ - flag, keep frequency ranges that are flagged.
+ - move over the raw data at full time resolution; FFT, replace with 0, mean or median; inverse FFT
+
+ or
+
+ - do not integrate, but move over raw data in full time resolution
+ - do fft
+ - flag on this data only
+ - replace with 0, mean or median
+ - inverse fft
+
+ In all these cases replacing with median would be best, but expensive.
+ Also, which median? compute it once on the raw data for all samples, or once per fft?
+
+ Option 1 is cheaper, since we flag only once, instead of integrationFactor times.
+ It may also be better due to the improved signal-to-noise ratio.
+*/
+
+namespace LOFAR {
+namespace RTCP {
+
+PreCorrelationNoChannelsFlagger::PreCorrelationNoChannelsFlagger(const Parset& parset, const unsigned nrStations, const unsigned nrSubbands, const unsigned nrChannels,
+ const unsigned nrSamplesPerIntegration, const float cutoffThreshold)
+:
+ Flagger(parset, nrStations, nrSubbands, nrChannels, cutoffThreshold, /*baseSentitivity*/ 1.0f,
+ getFlaggerStatisticsType(parset.onlinePreCorrelationFlaggingStatisticsType(getFlaggerStatisticsTypeString(FLAGGER_STATISTICS_WINSORIZED)))),
+ itsNrSamplesPerIntegration(nrSamplesPerIntegration)
+{
+ assert(itsNrSamplesPerIntegration % itsFFTSize == 0);
+ assert(nrChannels == 1);
+
+ itsIntegrationFactor = itsNrSamplesPerIntegration / itsFFTSize;
+
+ itsSamples.resize(itsFFTSize);
+ itsPowers.resize(itsFFTSize);
+ itsFlagsTime.resize(itsFFTSize);
+ itsFlagsFrequency.resize(itsFFTSize);
+ itsFFTBuffer.resize(itsFFTSize);
+}
+
+
+void PreCorrelationNoChannelsFlagger::initFFT()
+{
+#if defined HAVE_FFTW3
+ itsFFTWforwardPlan = fftwf_plan_dft_1d(itsFFTSize, (fftwf_complex *) &itsSamples[0], (fftwf_complex *) &itsFFTBuffer[0], FFTW_FORWARD, FFTW_MEASURE);
+ itsFFTWbackwardPlan = fftwf_plan_dft_1d(itsFFTSize, (fftwf_complex *) &itsFFTBuffer[0], (fftwf_complex *) &itsSamples[0], FFTW_FORWARD, FFTW_MEASURE);
+#elif defined HAVE_FFTW2
+#error TODO // @@@ TODO
+// itsFFTWforwardPlan = fftw_create_plan_specific(itsFFTsize, FFTW_FORWARD, FFTW_ESTIMATE | FFTW_USE_WISDOM, (fftw_complex *) &itsSamples[0], 2, (fftw_complex *) &itsFFTBuffer[0], 1);
+// itsFFTWbackwardPlan = fftw_create_plan_specific(itsFFTsize, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_USE_WISDOM, (fftw_complex *) &itsFFTBuffer[0], 1, (fftw_complex *) &itsSamples[0], 2);
+#endif
+}
+
+
+void PreCorrelationNoChannelsFlagger::forwardFFT()
+{
+#if defined HAVE_FFTW3
+ fftwf_execute(itsFFTWforwardPlan);
+#elif defined HAVE_FFTW2
+ fftw(itsFFTWforwardPlan, 2, (fftw_complex *) data, 2, 1, (fftw_complex *) &itsFFTedBuffer[0][0], 1, itsFFTsize);
+#endif
+}
+
+
+void PreCorrelationNoChannelsFlagger::backwardFFT()
+{
+#if defined HAVE_FFTW3
+ fftwf_execute(itsFFTWbackwardPlan);
+#elif defined HAVE_FFTW2
+ fftw(itsFFTWbackwardPlan, 2, (fftw_complex *) &itsFFTedBuffer[0][0], 1, itsFFTsize, (fftw_complex *) data, 2, 1);
+#endif
+}
+
+
+void PreCorrelationNoChannelsFlagger::flag(FilteredData* filteredData, unsigned currentSubband)
+{
+ NSTimer flaggerTimer("RFI noChannels flagger total", true, true);
+ NSTimer flaggerTimeTimer("RFI noChannels time flagger", true, true);
+ NSTimer flaggerFrequencyTimer("RFI noChannels time flagger", true, true);
+
+ flaggerTimer.start();
+
+ for(unsigned station = 0; station < itsNrStations; station++) {
+ initFlagsTime(station, filteredData); // copy flags to my local format
+ filteredData->resetFlags(); // Wipe original flags
+
+ // init frequency flags
+ for (unsigned i = 0; i < itsFFTSize; i++) {
+ itsFlagsFrequency[i] = false;
+ }
+
+ for (unsigned pol = 0; pol < NR_POLARIZATIONS; pol++) {
+ integrateAndCalculatePowers(station, pol, filteredData);
+
+ flaggerTimeTimer.start();
+ sumThresholdFlagger1D(itsPowers, itsFlagsTime, itsBaseSensitivity); // flag in time direction
+ flaggerTimeTimer.stop();
+
+ flaggerFrequencyTimer.start();
+ forwardFFT();
+
+ for (unsigned i = 0; i < itsFFTSize; i++) { // compute powers from FFT-ed data
+ fcomplex sample = itsFFTBuffer[i];
+ float power = real(sample) * real(sample) + imag(sample) * imag(sample);
+ itsPowers[i] = power;
+ }
+
+ sumThresholdFlagger1D(itsPowers, itsFlagsFrequency, itsBaseSensitivity); // flag in freq direction
+ flaggerFrequencyTimer.stop();
+ }
+
+ applyFlagsTime(station, filteredData); // copy flags from my original format into FilteredData again.
+ applyFlagsFrequency(station, filteredData); // do forward FFT; fix samples; backward FFT on the original samples in full resolution
+ }
+
+ flaggerTimer.stop();
+}
+
+
+void PreCorrelationNoChannelsFlagger::integrateAndCalculatePowers(unsigned station, unsigned pol, FilteredData* filteredData)
+{
+ for(unsigned i=0; i<itsFFTSize; i++) {
+ itsSamples[i] = makefcomplex(0, 0);
+ }
+
+ for(unsigned t=0; t<itsNrSamplesPerIntegration; t++) {
+ itsSamples[t/itsIntegrationFactor] += filteredData->samples[0][station][t][pol];
+ }
+
+ for (unsigned i = 0; i < itsFFTSize; i++) {
+ fcomplex sample = itsSamples[i];
+ float power = real(sample) * real(sample) + imag(sample) * imag(sample);
+ itsPowers[i] = power;
+ }
+}
+
+
+void PreCorrelationNoChannelsFlagger::initFlagsTime(unsigned station, FilteredData* filteredData)
+{
+ for (unsigned i = 0; i < itsFFTSize; i++) {
+ itsFlagsTime[i] = false;
+ }
+
+ // Use the original flags to initialize the flags.
+ // This could be done much faster by just iterating over the windows in the sparse flags set. // TODO
+ for (unsigned time = 0; time < itsNrSamplesPerIntegration; time++) {
+ if(filteredData->flags[0][station].test(time)) {
+ itsFlagsTime[time/itsIntegrationFactor] = true;
+ }
+ }
+}
+
+
+void PreCorrelationNoChannelsFlagger::applyFlagsTime(unsigned station, FilteredData* filteredData)
+{
+ const fcomplex zero = makefcomplex(0, 0);
+
+ for (unsigned i = 0; i < itsFFTSize; i++) {
+ if(itsFlagsTime[i]) {
+ unsigned startIndex = i * itsIntegrationFactor;
+
+ filteredData->flags[0][station].include(startIndex, startIndex+itsIntegrationFactor);
+
+ for (unsigned time = 0; time < itsIntegrationFactor; time++) {
+ unsigned globalIndex = i * itsIntegrationFactor + time;
+ for (unsigned pol = 0; pol < NR_POLARIZATIONS; pol++) {
+ filteredData->samples[0][station][globalIndex][pol] = zero;
+ }
+ }
+ }
+ }
+}
+
+// Do forward FFT; fix samples; backward FFT on the original samples in full resolution. Flags are already set in itsFlagsFrequency.
+// FFT followed by an inverse FFT multiplies all samples by N. Thus, we have to divide by N after we are done.
+void PreCorrelationNoChannelsFlagger::applyFlagsFrequency(unsigned station, FilteredData* filteredData)
+{
+ const fcomplex zero = makefcomplex(0, 0);
+
+ for (unsigned time = 0; time < itsIntegrationFactor; time++) {
+ unsigned startIndex = time * itsFFTSize;
+ for (unsigned pol = 0; pol < NR_POLARIZATIONS; pol++) {
+ for(unsigned minorTime=0; minorTime < itsFFTSize; minorTime++) {
+ itsSamples[minorTime] = filteredData->samples[0][station][startIndex+minorTime][pol];
+ }
+ forwardFFT();
+ for(unsigned minorTime=0; minorTime < itsFFTSize; minorTime++) {
+ if(itsFlagsFrequency[minorTime]) {
+ itsFFTBuffer[minorTime] = zero;
+ }
+ }
+ backwardFFT();
+ for(unsigned minorTime=0; minorTime < itsFFTSize; minorTime++) {
+ filteredData->samples[0][station][startIndex+minorTime][pol] = makefcomplex(real(itsSamples[minorTime]) / itsFFTSize, imag(itsSamples[minorTime]) / itsFFTSize);
+ }
+ }
+ }
+}
+
+fcomplex PreCorrelationNoChannelsFlagger::computeMedianSample(unsigned station, FilteredData* filteredData)
+{
+ // we have to copy the vector, nth_element changes the ordering, also, we want the median of both polarizations
+ std::vector<fcomplex> copy(itsNrSamplesPerIntegration * NR_POLARIZATIONS);
+ memcpy(copy.data(), &filteredData->samples[0][station][0][0], itsNrSamplesPerIntegration * NR_POLARIZATIONS * sizeof(fcomplex));
+
+ std::vector<float> powers(itsNrSamplesPerIntegration * NR_POLARIZATIONS);
+ for(unsigned i=0; i<itsNrSamplesPerIntegration * NR_POLARIZATIONS; i++) {
+ fcomplex sample = copy[i];
+ powers[i] = real(sample) * real(sample) + imag(sample) * imag(sample);
+ }
+
+ // calculate median, expensive, but nth_element is guaranteed to be O(n)
+ std::vector<float>::iterator it = powers.begin() + (powers.size() / 2);
+ std::nth_element(powers.begin(), it, powers.end());
+
+ float median = *it;
+
+ for(unsigned i=0; i<itsNrSamplesPerIntegration * NR_POLARIZATIONS; i++) {
+ if(powers[i] == median) {
+ return filteredData->samples[0][station][i/NR_POLARIZATIONS][i%NR_POLARIZATIONS];
+ }
+ }
+
+ return makefcomplex(0, 0);
+}
+
+
+PreCorrelationNoChannelsFlagger::~PreCorrelationNoChannelsFlagger()
+{
+#if defined HAVE_FFTW3
+ if(itsFFTWforwardPlan != 0) {
+ fftwf_destroy_plan(itsFFTWforwardPlan);
+ }
+ if(itsFFTWbackwardPlan != 0) {
+ fftwf_destroy_plan(itsFFTWbackwardPlan);
+ }
+#else
+ if(itsFFTWforwardPlan != 0) {
+ fftw_destroy_plan(itsFFTWforwardPlan);
+ }
+ if(itsFFTWbackwardPlan != 0) {
+ fftw_destroy_plan(itsFFTWbackwardPlan);
+ }
+#endif
+}
+
+} // namespace RTCP
+} // namespace LOFAR
diff --git a/RTCP/CNProc/src/PreCorrelationNoChannelsFlagger.h b/RTCP/CNProc/src/PreCorrelationNoChannelsFlagger.h
new file mode 100644
index 00000000000..b1014915863
--- /dev/null
+++ b/RTCP/CNProc/src/PreCorrelationNoChannelsFlagger.h
@@ -0,0 +1,63 @@
+#ifndef LOFAR_CNPROC_PRE_CORRELATION_NO_CHANNELS_FLAGGER_H
+#define LOFAR_CNPROC_PRE_CORRELATION_NO_CHANNELS_FLAGGER_H
+
+#include <Flagger.h>
+#include <Interface/FilteredData.h>
+
+#if defined HAVE_FFTW3
+#include <fftw3.h>
+#elif defined HAVE_FFTW2
+#include <fftw.h>
+#else
+#error Should have FFTW3 or FFTW2 installed
+#endif
+
+namespace LOFAR {
+namespace RTCP {
+
+// integrate in time untill we have itsFFTSize elements.
+// Flag on that in time direction.
+// Next, do FFT, flag in frequency direction, replace samples with median, inverseFFT
+class PreCorrelationNoChannelsFlagger : public Flagger {
+ public:
+ PreCorrelationNoChannelsFlagger(const Parset& parset, const unsigned nrStations, const unsigned nrSubbands, const unsigned nrChannels, const unsigned nrSamplesPerIntegration, float cutoffThreshold = 7.0f);
+
+ void flag(FilteredData* filteredData, unsigned currentSubband);
+
+ ~PreCorrelationNoChannelsFlagger();
+
+ private:
+
+ static const unsigned itsFFTSize = 256;
+
+ const unsigned itsNrSamplesPerIntegration;
+ unsigned itsIntegrationFactor;
+
+ void integrateAndCalculatePowers(unsigned station, unsigned pol, FilteredData* filteredData);
+ void initFlagsTime(unsigned station, FilteredData* filteredData);
+ void applyFlagsTime(unsigned station, FilteredData* filteredData);
+ void applyFlagsFrequency(unsigned station, FilteredData* filteredData);
+ fcomplex computeMedianSample(unsigned station, FilteredData* filteredData);
+
+ void initFFT();
+ void forwardFFT();
+ void backwardFFT();
+
+ vector<fcomplex> itsSamples; // [itsFFTSize]
+ vector<float> itsPowers; // [itsFFTSize]
+ vector<bool> itsFlagsTime; // [itsFFTSize]
+ vector<bool> itsFlagsFrequency; // [itsFFTSize]
+ vector<fcomplex> itsFFTBuffer; // [itsFFTSize]
+
+#if defined HAVE_FFTW3
+ fftwf_plan itsFFTWforwardPlan, itsFFTWbackwardPlan;
+#elif defined HAVE_FFTW2
+ fftw_plan itsFFTWforwardPlan, itsFFTWbackwardPlan;
+#endif
+
+};
+
+} // namespace RTCP
+} // namespace LOFAR
+
+#endif // LOFAR_CNPROC_PRE_CORRELATION_NO_CHANNELS_FLAGGER_H
diff --git a/RTCP/CNProc/test/FlaggerTest.parset b/RTCP/CNProc/test/FlaggerTest.parset
index a8aebb97fb7..9e5debbabf9 100644
--- a/RTCP/CNProc/test/FlaggerTest.parset
+++ b/RTCP/CNProc/test/FlaggerTest.parset
@@ -1,4 +1,4 @@
-OLAP.CNProc.integrationSteps = 12208
+OLAP.CNProc.integrationSteps = 12288
OLAP.CNProc.phaseOnePsets = [0..4]
OLAP.CNProc.phaseTwoPsets = [0..4]
OLAP.CNProc.phaseThreePsets = [0..4]
@@ -9,11 +9,13 @@ OLAP.CNProc.tabList = []
OLAP.CNProc.onlineFlagging = F # enable or disable all online flagging
-OLAP.CNProc.onlinePreCorrelationFlagging = T
+OLAP.CNProc.onlinePreCorrelationNoChannelsFlagging = F
+
+OLAP.CNProc.onlinePreCorrelationFlagging = F
#OLAP.CNProc.onlinePreCorrelationFlaggingType = INTEGRATED_THRESHOLD
OLAP.CNProc.onlinePreCorrelationFlaggingType = INTEGRATED_SUM_THRESHOLD_2D_WITH_HISTORY # THRESHOLD, INTEGRATED_THRESHOLD, INTEGRATED_THRESHOLD_2D, SUM_THRESHOLD, INTEGRATED_SUM_THRESHOLD, INTEGRATED_SUM_THRESHOLD_WITH_HISTORY, INTEGRATED_SMOOTHED_SUM_THRESHOLD, INTEGRATED_SMOOTHED_SUM_THRESHOLD_WITH_HISTORY, INTEGRATED_SUM_THRESHOLD_2D, INTEGRATED_SUM_THRESHOLD_2D_WITH_HISTORY
OLAP.CNProc.onlinePreCorrelationFlaggingStatisticsType = WINSORIZED # NORMAL, WINSORIZED
-OLAP.CNProc.onlinePostCorrelationFlaggingIntegration = 763
+OLAP.CNProc.onlinePostCorrelationFlaggingIntegration = 768
OLAP.CNProc.onlinePostCorrelationFlagging= F
OLAP.CNProc.onlinePostCorrelationFlaggingType = SMOOTHED_SUM_THRESHOLD_WITH_HISTORY # THRESHOLD, SUM_THRESHOLD, SMOOTHED_SUM_THRESHOLD, SMOOTHED_SUM_THRESHOLD_WITH_HISTORY
@@ -46,8 +48,8 @@ Observation.Beam[0].TiedArrayBeam[0].dispersionMeasure = 0
OLAP.IONProc.integrationSteps = 1
OLAP.CNProc_CoherentStokes.timeIntegrationFactor = 1
OLAP.CNProc_IncoherentStokes.timeIntegrationFactor = 1
-OLAP.CNProc_CoherentStokes.channelsPerSubband = 16
-OLAP.CNProc_IncoherentStokes.channelsPerSubband = 16
+OLAP.CNProc_CoherentStokes.channelsPerSubband = 32
+OLAP.CNProc_IncoherentStokes.channelsPerSubband = 32
#OLAP.CNProc_CoherentStokes.which = IQUV
#OLAP.CNProc_IncoherentStokes.which = IQUV
OLAP.CNProc_CoherentStokes.which = I
@@ -74,7 +76,7 @@ Observation.subbandList = [200..231]
Observation.beamList = [32*0]
Observation.rspBoardList = [32*0]
Observation.rspSlotList = [0..31]
-Observation.channelsPerSubband = 16
+Observation.channelsPerSubband = 32
Observation.sampleClock = 200
Observation.nrSlotsInFrame = 32
Observation.ObsID = 1000000
diff --git a/RTCP/CNProc/test/RFI-test.parset b/RTCP/CNProc/test/RFI-test.parset
index aeecbe565e0..3c5a01111d7 100644
--- a/RTCP/CNProc/test/RFI-test.parset
+++ b/RTCP/CNProc/test/RFI-test.parset
@@ -9,7 +9,9 @@ OLAP.CNProc.tabList = []
OLAP.CNProc.onlineFlagging = T # enable or disable all online flagging
-OLAP.CNProc.onlinePreCorrelationFlagging = T
+OLAP.CNProc.onlinePreCorrelationNoChannelsFlagging = T
+
+OLAP.CNProc.onlinePreCorrelationFlagging = F
OLAP.CNProc.onlinePreCorrelationFlaggingType = INTEGRATED_SMOOTHED_SUM_THRESHOLD_WITH_HISTORY # THRESHOLD, INTEGRATED_THRESHOLD, SUM_THRESHOLD, INTEGRATED_SUM_THRESHOLD, INTEGRATED_SMOOTHED_SUM_THRESHOLD, INTEGRATED_SMOOTHED_SUM_THRESHOLD_WITH_HISTORY
OLAP.CNProc.onlinePreCorrelationFlaggingStatisticsType = WINSORIZED # NORMAL, WINSORIZED
diff --git a/RTCP/CNProc/test/tInversePPF.cc b/RTCP/CNProc/test/tInversePPF.cc
index 27c8b321edf..1217b19d129 100644
--- a/RTCP/CNProc/test/tInversePPF.cc
+++ b/RTCP/CNProc/test/tInversePPF.cc
@@ -52,7 +52,8 @@ const static unsigned nrTaps = 16;
static unsigned nrSubbands = 248;
//static unsigned nrSubbands = 4;
static unsigned nrChannels = 1; // for the NuMoon pipeline, there are no separate channels.
-static unsigned nrSamplesPerIntegration = 768 * 256 / 4; // one quarter of a second
+//static unsigned nrSamplesPerIntegration = 768 * 256 / 4; // one quarter of a second
+static unsigned nrSamplesPerIntegration = 19648; // roughly 0.1 seconds
//static unsigned nrSamplesPerIntegration = 64;
static double sampleRate = 195312.5;
static double centerFrequency = (nrSamplesPerIntegration / 2) * sampleRate;
@@ -110,7 +111,7 @@ static void performStationFFT(TransposedBeamFormedData& transposedBeamFormedData
// Put data in the right order, go from half complex to normal format
for (unsigned subbandIndex = 0; subbandIndex < subbandList.size(); subbandIndex++) {
unsigned subband = subbandList[subbandIndex];
- fcomplex sample = makefcomplex(fftOutData[subband], fftOutData[onStationFilterSize - subband]);
+ fcomplex sample = makefcomplex(fftOutData[subband], fftOutData[onStationFilterSize - subband - 1]);
transposedBeamFormedData.samples[subband][0 /* channel, but there is only one now */][time] = sample;
}
}
@@ -180,7 +181,7 @@ static void fftTest() {
#if 0
// Put data in the right order, go from half complex to normal format
for (unsigned subband = 0; subband < nrSubbands; subband++) {
- fcomplex sample = makefcomplex(fftOutData[subband], fftOutData[onStationFilterSize - subband]);
+ fcomplex sample = makefcomplex(fftOutData[subband], fftOutData[onStationFilterSize - subband - 1]);
transposedBeamFormedData.samples[subband][0 /* channel */][time] = sample;
}
#endif
@@ -190,12 +191,10 @@ static void fftTest() {
float maxError = 0.0f;
for (unsigned time = 0; time < onStationFilterSize; time++) {
-
float error = fabsf(inputData[time] - (fftInData[time]/((float)onStationFilterSize))); // the error
if(error > maxError) {
maxError = error;
}
-
// fprintf(stdout, "%20.10lf\n", error);
// fprintf(stdout, "%20.10lf\n", fftInData[time]);
}
diff --git a/RTCP/Interface/include/Interface/Parset.h b/RTCP/Interface/include/Interface/Parset.h
index 86117d77816..8b049dd76c7 100644
--- a/RTCP/Interface/include/Interface/Parset.h
+++ b/RTCP/Interface/include/Interface/Parset.h
@@ -167,6 +167,7 @@ class Parset: public ParameterSet
bool onlineFlagging() const;
bool onlinePreCorrelationFlagging() const;
+ bool onlinePreCorrelationNoChannelsFlagging() const;
bool onlinePostCorrelationFlagging() const;
bool onlinePostCorrelationFlaggingDetectBrokenStations() const;
unsigned onlinePreCorrelationFlaggingIntegration() const;
@@ -780,6 +781,11 @@ inline bool Parset::onlinePreCorrelationFlagging() const
return getBool("OLAP.CNProc.onlinePreCorrelationFlagging", false);
}
+inline bool Parset::onlinePreCorrelationNoChannelsFlagging() const
+{
+ return getBool("OLAP.CNProc.onlinePreCorrelationNoChannelsFlagging", false);
+}
+
inline bool Parset::onlinePostCorrelationFlagging() const
{
return getBool("OLAP.CNProc.onlinePostCorrelationFlagging", false);
--
GitLab