diff --git a/.gitattributes b/.gitattributes index f420ec09b7a6cad9ed5be39382a175c6cf6f15ff..0236be50304cef6b8922c17d389090efaa3a146c 100644 --- a/.gitattributes +++ b/.gitattributes @@ -1991,6 +1991,10 @@ RTCP/CNProc/src/BeamFormerAsm3St6Bm.inc -text RTCP/CNProc/src/BeamFormerAsm6St3Bm.inc -text RTCP/CNProc/src/CN_Math.h -text RTCP/CNProc/src/CN_Processing.machinefile -text +RTCP/CNProc/src/Dedispersion.cc -text +RTCP/CNProc/src/Dedispersion.h -text +RTCP/CNProc/src/DedispersionAsm.S -text +RTCP/CNProc/src/DedispersionAsm.h -text RTCP/CNProc/src/FIR_InvertedStationPPFWeights.h -text RTCP/CNProc/src/FIR_OriginalCepPPFWeights.h -text RTCP/CNProc/src/FIR_OriginalStationPPFWeights.h -text @@ -2013,6 +2017,7 @@ RTCP/CNProc/test/filterTestOutput.dat -text RTCP/CNProc/test/filterTestResult.ps -text RTCP/CNProc/test/inversePPFTestOutput.dat -text RTCP/CNProc/test/inversePPFTestResult.ps -text +RTCP/CNProc/test/tDedispersion.cc -text RTCP/CNProc/test/tInversePPF.cc -text RTCP/CNProc/test/tPencilBeamFormer.cc -text RTCP/CNProc/test/tPencilBeamFormer.sh -text diff --git a/RTCP/CNProc/src/CMakeLists.txt b/RTCP/CNProc/src/CMakeLists.txt index 2e08f8ee239dd98651580aead118d7db040ece30..1dcdbfdc2da6cdd314f9a3165ed5a8b9efbace88 100644 --- a/RTCP/CNProc/src/CMakeLists.txt +++ b/RTCP/CNProc/src/CMakeLists.txt @@ -20,10 +20,11 @@ set(cnproc_LIB_SRCS BeamFormer.cc CN_Processing.cc Correlator.cc + Dedispersion.cc FCNP_ClientStream.cc FIR.cc - InversePPF.cc FilterBank.cc + InversePPF.cc LocationInfo.cc PPF.cc Flagger.cc @@ -38,6 +39,7 @@ if(CMAKE_ASM-BGP_COMPILER_WORKS) list(APPEND cnproc_LIB_SRCS BeamFormerAsm.S CorrelatorAsm.S + DedispersionAsm.S FIR_Asm.S FFT_Asm.S StokesAsm.S) diff --git a/RTCP/CNProc/src/CN_Processing.cc b/RTCP/CNProc/src/CN_Processing.cc index 3de9075bb27a8d3b8a719eafae70ae4bb7b16400..b151ae79394fec163039bf006b2520cc01a5d745 100644 --- a/RTCP/CNProc/src/CN_Processing.cc +++ b/RTCP/CNProc/src/CN_Processing.cc @@ -94,6 +94,8 @@ template <typename SAMPLE_TYPE> CN_Processing<SAMPLE_TYPE>::CN_Processing(Stream itsCoherentStokes(0), itsIncoherentStokes(0), itsCorrelator(0), + itsDedispersionBeforeBeamForming(0), + itsDedispersionAfterBeamForming(0), itsDoOnlineFlagging(0), itsPreCorrelationFlagger(0), itsPostCorrelationFlagger(0) @@ -255,17 +257,21 @@ template <typename SAMPLE_TYPE> void CN_Processing<SAMPLE_TYPE>::preprocess(CN_C itsPPF = new PPF<SAMPLE_TYPE>(itsNrStations, nrChannels, nrSamplesPerIntegration, configuration.sampleRate() / nrChannels, configuration.delayCompensation(), configuration.correctBandPass(), itsLocationInfo.rank() == 0); - if(itsDoOnlineFlagging) { + if (configuration.dispersionMeasure() != 0) + if (configuration.outputIncoherentStokes() || configuration.outputCorrelatedData() || itsNrBeamFormedStations < itsNrBeams) + itsDedispersionBeforeBeamForming = new DedispersionBeforeBeamForming(configuration, itsPlan->itsFilteredData, itsCurrentSubband->list()); + else + itsDedispersionAfterBeamForming = new DedispersionAfterBeamForming(configuration, itsPlan->itsBeamFormedData, itsCurrentSubband->list()); + + if(itsDoOnlineFlagging) itsPreCorrelationFlagger = new PreCorrelationFlagger(itsNrStations, nrChannels, nrSamplesPerIntegration); - } itsIncoherentStokes = new Stokes(itsNrStokes, nrChannels, nrSamplesPerIntegration, nrSamplesPerStokesIntegration, configuration.stokesNrChannelsPerSubband() ); itsCorrelator = new Correlator(itsBeamFormer->getStationMapping(), nrChannels, nrSamplesPerIntegration); - if(itsDoOnlineFlagging) { + if (itsDoOnlineFlagging) itsPostCorrelationFlagger = new PostCorrelationFlagger(itsNrStations, nrChannels); - } } if (itsHasPhaseThree && itsPhaseThreeDisjunct) { @@ -495,6 +501,7 @@ template <typename SAMPLE_TYPE> int CN_Processing<SAMPLE_TYPE>::transposeBeams(u return beamToProcess ? myBeam : -1; } + template <typename SAMPLE_TYPE> void CN_Processing<SAMPLE_TYPE>::filter() { #if defined HAVE_MPI @@ -524,11 +531,51 @@ template <typename SAMPLE_TYPE> void CN_Processing<SAMPLE_TYPE>::filter() #endif // HAVE_MPI } + +template <typename SAMPLE_TYPE> void CN_Processing<SAMPLE_TYPE>::dedisperseBeforeBeamForming() +{ + if (itsDedispersionBeforeBeamForming != 0) { +#if defined HAVE_MPI + if (LOG_CONDITION) + LOG_DEBUG_STR(itsLogPrefix << "Start dedispersion at " << MPI_Wtime()); +#endif + + static NSTimer timer("dedispersion (before BF) timer", true, true); + + computeTimer.start(); + timer.start(); + itsDedispersionBeforeBeamForming->dedisperse(itsPlan->itsFilteredData, *itsCurrentSubband); + timer.stop(); + computeTimer.stop(); + } +} + + +template <typename SAMPLE_TYPE> void CN_Processing<SAMPLE_TYPE>::dedisperseAfterBeamForming() +{ + if (itsDedispersionAfterBeamForming != 0) { +#if defined HAVE_MPI + if (LOG_CONDITION) + LOG_DEBUG_STR(itsLogPrefix << "Start dedispersion at " << MPI_Wtime()); +#endif + + static NSTimer timer("dedispersion (after BF) timer", true, true); + + computeTimer.start(); + timer.start(); + itsDedispersionAfterBeamForming->dedisperse(itsPlan->itsBeamFormedData, *itsCurrentSubband); + timer.stop(); + computeTimer.stop(); + } +} + + template <typename SAMPLE_TYPE> void CN_Processing<SAMPLE_TYPE>::preCorrelationFlagging() { itsPreCorrelationFlagger->flag(itsPlan->itsFilteredData); } + template <typename SAMPLE_TYPE> void CN_Processing<SAMPLE_TYPE>::mergeStations() { #if defined HAVE_MPI @@ -788,6 +835,7 @@ template <typename SAMPLE_TYPE> void CN_Processing<SAMPLE_TYPE>::process(unsigne if( itsPlan->calculate( itsPlan->itsFilteredData ) ) { filter(); mergeStations(); // create superstations + dedisperseBeforeBeamForming(); } if(itsDoOnlineFlagging) { @@ -796,6 +844,7 @@ template <typename SAMPLE_TYPE> void CN_Processing<SAMPLE_TYPE>::process(unsigne if( itsPlan->calculate( itsPlan->itsBeamFormedData ) ) { formBeams(); + dedisperseAfterBeamForming(); } if( itsPlan->calculate( itsPlan->itsCorrelatedData ) ) { @@ -869,23 +918,25 @@ template <typename SAMPLE_TYPE> void CN_Processing<SAMPLE_TYPE>::postprocess() { #if defined HAVE_MPI - delete itsAsyncTransposeInput; itsAsyncTransposeInput = 0; - delete itsAsyncTransposeBeams; itsAsyncTransposeBeams = 0; + delete itsAsyncTransposeInput; itsAsyncTransposeInput = 0; + delete itsAsyncTransposeBeams; itsAsyncTransposeBeams = 0; #endif // HAVE_MPI - delete itsBeamFormer; itsBeamFormer = 0; - delete itsPPF; itsPPF = 0; - delete itsPreCorrelationFlagger; itsPreCorrelationFlagger = 0; - delete itsPostCorrelationFlagger; itsPostCorrelationFlagger = 0; - delete itsCorrelator; itsCorrelator = 0; - delete itsCoherentStokes; itsCoherentStokes = 0; - delete itsIncoherentStokes; itsIncoherentStokes = 0; - - delete itsCurrentSubband; itsCurrentSubband = 0; - delete itsCurrentBeam; itsCurrentBeam = 0; - - for (unsigned i = 0; i < itsOutputStreams.size(); i++) { + delete itsBeamFormer; itsBeamFormer = 0; + delete itsPPF; itsPPF = 0; + delete itsDedispersionBeforeBeamForming; itsDedispersionBeforeBeamForming = 0; + delete itsDedispersionAfterBeamForming; itsDedispersionAfterBeamForming = 0; + delete itsCorrelator; itsCorrelator = 0; + delete itsCoherentStokes; itsCoherentStokes = 0; + delete itsIncoherentStokes; itsIncoherentStokes = 0; + delete itsPreCorrelationFlagger; itsPreCorrelationFlagger = 0; + delete itsPostCorrelationFlagger; itsPostCorrelationFlagger = 0; + + delete itsCurrentSubband; itsCurrentSubband = 0; + delete itsCurrentBeam; itsCurrentBeam = 0; + + for (unsigned i = 0; i < itsOutputStreams.size(); i++) delete itsOutputStreams[i]; - } + itsOutputStreams.clear(); delete itsPlan; itsPlan = 0; diff --git a/RTCP/CNProc/src/CN_Processing.h b/RTCP/CNProc/src/CN_Processing.h index 823a8359ad032842d66039eb3ed27591eab345ea..08c0066f81ca0ff27506a7b5cfa4f306a9f44c95 100644 --- a/RTCP/CNProc/src/CN_Processing.h +++ b/RTCP/CNProc/src/CN_Processing.h @@ -40,6 +40,7 @@ #include <AsyncTranspose.h> #include <AsyncTransposeBeams.h> #include <BeamFormer.h> +#include <Dedispersion.h> #include <PPF.h> #include <Correlator.h> #include <Stokes.h> @@ -81,6 +82,8 @@ template <typename SAMPLE_TYPE> class CN_Processing : public CN_Processing_Base, void transposeInput(); int transposeBeams(unsigned block); void filter(); + void dedisperseBeforeBeamForming(); + void dedisperseAfterBeamForming(); void preCorrelationFlagging(); void mergeStations(); void formBeams(); @@ -142,6 +145,8 @@ template <typename SAMPLE_TYPE> class CN_Processing : public CN_Processing_Base, BeamFormer *itsBeamFormer; Stokes *itsCoherentStokes, *itsIncoherentStokes; Correlator *itsCorrelator; + DedispersionBeforeBeamForming *itsDedispersionBeforeBeamForming; + DedispersionAfterBeamForming *itsDedispersionAfterBeamForming; bool itsDoOnlineFlagging; PreCorrelationFlagger *itsPreCorrelationFlagger; PostCorrelationFlagger *itsPostCorrelationFlagger; diff --git a/RTCP/CNProc/src/Dedispersion.cc b/RTCP/CNProc/src/Dedispersion.cc new file mode 100644 index 0000000000000000000000000000000000000000..f7c5fb0291b1c94420c430f0f427a1bd897dc195 --- /dev/null +++ b/RTCP/CNProc/src/Dedispersion.cc @@ -0,0 +1,193 @@ +//# Always #include <lofar_config.h> first! +#include <lofar_config.h> + +#include <CN_Math.h> +#include <Dedispersion.h> +#include <DedispersionAsm.h> +#include <Common/Timer.h> + +#include <algorithm> + + +#if defined HAVE_FFTW3 +#include <fftw3.h> +#include <vector> +#elif defined HAVE_FFTW2 +#include <fftw.h> +#else +#error Should have FFTW3 or FFTW2 installed +#endif + + +namespace LOFAR { +namespace RTCP { + + +Dedispersion::Dedispersion(CN_Configuration &configuration, const std::vector<unsigned> &subbands) +: + itsNrChannels(configuration.nrChannelsPerSubband()), + itsNrSamplesPerIntegration(configuration.nrSamplesPerIntegration()), + itsFFTsize(configuration.dedispersionFFTsize()), + itsChannelBandwidth(configuration.sampleRate() / itsNrChannels), + itsFFTedBuffer(NR_POLARIZATIONS, itsFFTsize) +{ + initChirp(configuration, subbands); +} + + +DedispersionBeforeBeamForming::DedispersionBeforeBeamForming(CN_Configuration &configuration, FilteredData *filteredData, const std::vector<unsigned> &subbands) +: + Dedispersion(configuration, subbands), + itsNrStations(configuration.nrMergedStations()) +{ + initFFT(&filteredData->samples[0][0][0][0]); +} + + +DedispersionAfterBeamForming::DedispersionAfterBeamForming(CN_Configuration &configuration, BeamFormedData *beamFormedData, const std::vector<unsigned> &subbands) +: + Dedispersion(configuration, subbands), + itsNrBeams(configuration.flysEye() ? configuration.nrMergedStations() : configuration.nrPencilBeams()) +{ + initFFT(&beamFormedData->samples[0][0][0][0]); +} + + +Dedispersion::~Dedispersion() +{ +#if defined HAVE_FFTW3 + fftwf_destroy_plan(itsFFTWforwardPlan); + fftwf_destroy_plan(itsFFTWbackwardPlan); +#else + fftw_destroy_plan(itsFFTWforwardPlan); + fftw_destroy_plan(itsFFTWbackwardPlan); +#endif + + for (unsigned i = 0; i < itsChirp.size(); i ++) + delete itsChirp[i]; +} + + +void Dedispersion::initFFT(fcomplex *data) +{ +#if defined HAVE_FFTW3 + itsFFTWforwardPlan = fftwf_plan_many_dft(1, (int *) &itsFFTsize, 2, (fftwf_complex *) data, 0, 2, 1, (fftwf_complex *) &itsFFTedBuffer[0][0], 0, 1, itsFFTsize, FFTW_FORWARD, FFTW_MEASURE); + itsFFTWbackwardPlan = fftwf_plan_many_dft(1, (int *) &itsFFTsize, 2, (fftwf_complex *) &itsFFTedBuffer[0][0], 0, 1, itsFFTsize, (fftwf_complex *) data, 0, 2, 1, FFTW_BACKWARD, FFTW_MEASURE); +#elif defined HAVE_FFTW2 +#if defined HAVE_BGP + fftw_import_wisdom_from_string("(FFTW-2.1.5 (196608 529 1 0 1 2 1 77 0) (98304 529 1 0 1 2 1 99 0) (49152 529 1 0 1 2 1 715 0) (24576 529 1 0 1 2 1 715 0) (12288 529 1 0 1 2 1 715 0) (6144 529 1 0 1 2 1 77 0) (3072 529 1 0 1 2 1 715 0) (1536 529 1 0 1 2 1 187 0) (768 529 1 0 1 2 1 143 0) (384 529 1 0 1 2 1 143 0) (192 529 1 0 1 2 1 143 0) (96 529 1 0 1 2 1 143 0) (48 529 1 0 1 2 1 143 0) (24 529 1 0 1 2 1 143 0) (12 529 1 0 1 2 0 276 0) (6 529 1 0 1 2 0 144 0) (3 529 1 0 1 2 0 78 0) (196608 529 -1 0 2 1 1 704 0) (98304 529 -1 0 2 1 1 704 0) (49152 529 -1 0 2 1 1 704 0) (24576 529 -1 0 2 1 1 704 0) (12288 529 -1 0 2 1 1 704 0) (6144 529 -1 0 2 1 1 704 0) (3072 529 -1 0 2 1 1 132 0) (1536 529 -1 0 2 1 1 132 0) (768 529 -1 0 2 1 1 132 0) (384 529 -1 0 2 1 1 132 0) (192 529 -1 0 2 1 1 352 0) (96 529 -1 0 2 1 1 132 0) (48 529 -1 0 2 1 1 132 0) (24 529 -1 0 2 1 1 132 0) (12 529 -1 0 2 1 0 265 0) (6 529 -1 0 2 1 0 133 0) (3 529 -1 0 2 1 0 67 0) (2 529 -1 0 2 1 0 45 0) (4 529 -1 0 2 1 0 89 0) (8 529 -1 0 2 1 0 177 0) (16 529 -1 0 2 1 0 353 0) (32 529 -1 0 2 1 0 705 0) (64 529 -1 0 2 1 0 1409 0) (128 529 -1 0 2 1 0 2817 0) (256 529 -1 0 2 1 1 352 0) (512 529 -1 0 2 1 1 352 0) (1024 529 -1 0 2 1 1 704 0) (2048 529 -1 0 2 1 1 704 0) (4096 529 -1 0 2 1 1 704 0) (8192 529 -1 0 2 1 1 352 0) (16384 529 -1 0 2 1 1 704 0) (32768 529 -1 0 2 1 1 704 0) (65536 529 -1 0 2 1 1 704 0) (2 529 1 0 1 2 0 56 0) (4 529 1 0 1 2 0 100 0) (8 529 1 0 1 2 0 188 0) (16 529 1 0 1 2 0 364 0) (32 529 1 0 1 2 0 716 0) (64 529 1 0 1 2 0 1420 0) (128 529 1 0 1 2 0 2828 0) (256 529 1 0 1 2 1 715 0) (512 529 1 0 1 2 1 187 0) (1024 529 1 0 1 2 1 715 0) (2048 529 1 0 1 2 1 715 0) (4096 529 1 0 1 2 1 715 0) (8192 529 1 0 1 2 1 1419 0) (16384 529 1 0 1 2 1 99 0) (32768 529 1 0 1 2 1 715 0) (65536 529 1 0 1 2 1 715 0))"); +#endif + + itsFFTWforwardPlan = fftw_create_plan_specific(itsFFTsize, FFTW_FORWARD, FFTW_ESTIMATE | FFTW_USE_WISDOM, (fftw_complex *) data, 2, (fftw_complex *) &itsFFTedBuffer[0][0], 1); + itsFFTWbackwardPlan = fftw_create_plan_specific(itsFFTsize, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_USE_WISDOM, (fftw_complex *) &itsFFTedBuffer[0][0], 1, (fftw_complex *) data, 2); +#endif +} + + +void Dedispersion::forwardFFT(fcomplex *data) +{ +#if defined HAVE_FFTW3 + fftwf_execute_dft(itsFFTWforwardPlan, (fftwf_complex *) data, (fftwf_complex *) &itsFFTedBuffer[0][0]); +#elif defined HAVE_FFTW2 + fftw(itsFFTWforwardPlan, 2, (fftw_complex *) data, 2, 1, (fftw_complex *) &itsFFTedBuffer[0][0], 1, itsFFTsize); +#endif +} + + +void Dedispersion::backwardFFT(fcomplex *data) +{ +#if defined HAVE_FFTW3 + fftwf_execute_dft(itsFFTWbackwardPlan, (fftwf_complex *) &itsFFTedBuffer[0][0], (fftwf_complex *) data); +#elif defined HAVE_FFTW2 + fftw(itsFFTWbackwardPlan, 2, (fftw_complex *) &itsFFTedBuffer[0][0], 1, itsFFTsize, (fftw_complex *) data, 2, 1); +#endif +} + + +void Dedispersion::initChirp(CN_Configuration &configuration, const std::vector<unsigned> &subbands) +{ + itsChirp.resize(*std::max_element(subbands.begin(), subbands.end()) + 1, 0); +//std::cout << "newcurve linetype solid linethickness 3 marktype none color 0 .7 0 pts" << std::endl; + + for (unsigned subbandIndex = 0; subbandIndex < subbands.size(); subbandIndex ++) { + unsigned subband = subbands[subbandIndex]; + double subbandFrequency = configuration.refFreqs()[subbandIndex]; + double channel0frequency = subbandFrequency - (itsNrChannels * 0.5) * itsChannelBandwidth; + double binWidth = itsChannelBandwidth / itsFFTsize; + double dmConst = configuration.dispersionMeasure() * 2 * M_PI / 2.41e-16; + + itsChirp[subband] = new Matrix<fcomplex>(itsNrChannels, itsFFTsize); + + for (unsigned channel = 0; channel < itsNrChannels; channel ++) { + double channelFrequency = channel0frequency + channel * itsChannelBandwidth; + + for (unsigned n = 0; n < itsFFTsize; n ++) { + double binFrequency = n * binWidth; + + if (n > itsFFTsize / 2) + binFrequency -= itsChannelBandwidth; + + double frequencyDiv = binFrequency / channelFrequency; + double frequencyFac = frequencyDiv * frequencyDiv / (channelFrequency + binFrequency); + dcomplex dfactor = cosisin(dmConst * frequencyFac); + fcomplex factor = makefcomplex(real(dfactor), -imag(dfactor)); + float taper = sqrt(1 + pow(binFrequency / (.47 * itsChannelBandwidth), 80)); +//if (channel == 0) std::cout << n << ' ' << 1/taper << std::endl; + + (*itsChirp[subband])[channel][n] = factor / (taper * itsFFTsize); + } + } + } +} + + +void Dedispersion::applyChirp(unsigned subband, unsigned channel) +{ + static NSTimer chirpTimer("chirp timer", true, true); + const fcomplex *chirp = &(*itsChirp[subband])[channel][0]; + + chirpTimer.start(); + +#if defined HAVE_BGP + _apply_chirp(&itsFFTedBuffer[0][0], &itsFFTedBuffer[1][0], chirp, itsFFTsize); +#else + for (unsigned time = 0; time < itsFFTsize; time ++) { + itsFFTedBuffer[0][time] *= chirp[time]; + itsFFTedBuffer[1][time] *= chirp[time]; + } +#endif + + chirpTimer.stop(); +} + + +void DedispersionBeforeBeamForming::dedisperse(FilteredData *filteredData, unsigned subband) +{ + for (unsigned channel = 0; channel < itsNrChannels; channel ++) { + for (unsigned station = 0; station < itsNrStations; station ++) { + for (unsigned block = 0; block < itsNrSamplesPerIntegration; block += itsFFTsize) { + forwardFFT(&filteredData->samples[channel][station][block][0]); + applyChirp(subband, channel); + backwardFFT(&filteredData->samples[channel][station][block][0]); + } + } + } +} + + +void DedispersionAfterBeamForming::dedisperse(BeamFormedData *beamFormedData, unsigned subband) +{ + for (unsigned beam = 0; beam < itsNrBeams; beam ++) { + for (unsigned channel = 0; channel < itsNrChannels; channel ++) { + for (unsigned block = 0; block < itsNrSamplesPerIntegration; block += itsFFTsize) { + forwardFFT(&beamFormedData->samples[beam][channel][block][0]); + applyChirp(subband, channel); + backwardFFT(&beamFormedData->samples[beam][channel][block][0]); + } + } + } +} + + +} // namespace RTCP +} // namespace LOFAR diff --git a/RTCP/CNProc/src/Dedispersion.h b/RTCP/CNProc/src/Dedispersion.h new file mode 100644 index 0000000000000000000000000000000000000000..9d9495e4a8d2a6d0972de942ac11f645269d0414 --- /dev/null +++ b/RTCP/CNProc/src/Dedispersion.h @@ -0,0 +1,88 @@ +#ifndef LOFAR_CNPROC_DEDISPERSION_H +#define LOFAR_CNPROC_DEDISPERSION_H + + +#include <Common/lofar_complex.h> +#include <Interface/BeamFormedData.h> +#include <Interface/CN_Configuration.h> +#include <Interface/FilteredData.h> +#include <Interface/MultiDimArray.h> + +#include <vector> + + +#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 { + + +class Dedispersion +{ + protected: + Dedispersion(CN_Configuration &, const std::vector<unsigned> &subbands); + + public: + ~Dedispersion(); + + protected: + void initFFT(fcomplex *data); + void forwardFFT(fcomplex *data); + void backwardFFT(fcomplex *data); + + const unsigned itsNrChannels, itsNrSamplesPerIntegration, itsFFTsize; + const double itsChannelBandwidth; + + Matrix<fcomplex> itsFFTedBuffer; + +#if defined HAVE_FFTW3 + fftwf_plan itsFFTWforwardPlan, itsFFTWbackwardPlan; +#elif defined HAVE_FFTW2 + fftw_plan itsFFTWforwardPlan, itsFFTWbackwardPlan; +#endif + + void initChirp(CN_Configuration &configuration, const std::vector<unsigned> &subbands); + void applyChirp(unsigned subband, unsigned channel); + + std::vector<Matrix<fcomplex> *> itsChirp; // (*[subband])[channel][time] +}; + + +class DedispersionBeforeBeamForming : public Dedispersion +{ + public: + DedispersionBeforeBeamForming(CN_Configuration &, FilteredData *, const std::vector<unsigned> &subbands); + + void dedisperse(FilteredData *, unsigned subband); + + private: + const unsigned itsNrStations; +}; + + +class DedispersionAfterBeamForming : public Dedispersion +{ + public: + DedispersionAfterBeamForming(CN_Configuration &, BeamFormedData *, const std::vector<unsigned> &subbands); + + void dedisperse(BeamFormedData *, unsigned subband); + + private: + const unsigned itsNrBeams; +}; + + + + +} // namespace RTCP +} // namespace LOFAR + +#endif diff --git a/RTCP/CNProc/src/DedispersionAsm.S b/RTCP/CNProc/src/DedispersionAsm.S new file mode 100644 index 0000000000000000000000000000000000000000..ef61366229a7aa4dec023e7ad9153bb46510fade --- /dev/null +++ b/RTCP/CNProc/src/DedispersionAsm.S @@ -0,0 +1,122 @@ +#if defined HAVE_BGP + +.global _apply_chirp +_apply_chirp: + + li 12,-16 # push call-saved registers + stfpdux 14,1,12 + stfpdux 15,1,12 + stfpdux 16,1,12 + stfpdux 17,1,12 + stfpdux 18,1,12 + stfpdux 19,1,12 + + srwi 6,6,2 + subi 6,6,1 + mtctr 6 + li 8,8 + mr 9,3 + subi 10,4,8 + + lfpsx 0,0,3 + lfpsx 4,0,4 + lfpsx 8,0,5 + + lfpsux 1,3,8 + lfpsux 5,4,8 + lfpsux 9,5,8 + + lfpsux 2,3,8 + fxpmul 12,0,8 + lfpsux 6,4,8 + fxpmul 16,4,8 + lfpsux 10,5,8 + + lfpsux 3,3,8 + fxpmul 13,1,9 + lfpsux 7,4,8 + fxpmul 17,5,9 + lfpsux 11,5,8 + fxcxnpma 12,0,8,12 + fxcxnpma 16,4,8,16 + + lfpsux 0,3,8 + fxpmul 14,2,10 + lfpsux 4,4,8 + fxpmul 18,6,10 + lfpsux 8,5,8 + fxcxnpma 13,1,9,13 + stfpsx 12,0,9 + fxcxnpma 17,5,9,17 + +0: + lfpsux 1,3,8 + stfpsux 16,10,8 + fxpmul 15,3,11 + lfpsux 5,4,8 + fxpmul 19,7,11 + lfpsux 9,5,8 + fxcxnpma 14,2,10,14 + stfpsux 13,9,8 + fxcxnpma 18,6,10,18 + + lfpsux 2,3,8 + stfpsux 17,10,8 + fxpmul 12,0,8 + lfpsux 6,4,8 + fxpmul 16,4,8 + lfpsux 10,5,8 + fxcxnpma 15,3,11,15 + stfpsux 14,9,8 + fxcxnpma 19,7,11,19 + + lfpsux 3,3,8 + stfpsux 18,10,8 + fxpmul 13,1,9 + lfpsux 7,4,8 + fxpmul 17,5,9 + lfpsux 11,5,8 + fxcxnpma 12,0,8,12 + stfpsux 15,9,8 + fxcxnpma 16,4,8,16 + + lfpsux 0,3,8 + stfpsux 19,10,8 + fxpmul 14,2,10 + lfpsux 4,4,8 + fxpmul 18,6,10 + lfpsux 8,5,8 + fxcxnpma 13,1,9,13 + stfpsux 12,9,8 + fxcxnpma 17,5,9,17 + + bdnz 0b + + stfpsux 16,10,8 + fxpmul 15,3,11 + fxpmul 19,7,11 + fxcxnpma 14,2,10,14 + stfpsux 13,9,8 + fxcxnpma 18,6,10,18 + stfpsux 17,10,8 + + fxcxnpma 15,3,11,15 + stfpsux 14,9,8 + fxcxnpma 19,7,11,19 + stfpsux 18,10,8 + + stfpsux 15,9,8 + stfpsux 19,10,8 + + li 12,16 # restore call-saved registers + lfpdx 19,0,1 + lfpdux 18,1,12 + lfpdux 17,1,12 + lfpdux 16,1,12 + lfpdux 15,1,12 + lfpdux 14,1,12 + + addi 1,1,16 # reset stack pointer + + blr +#endif diff --git a/RTCP/CNProc/src/DedispersionAsm.h b/RTCP/CNProc/src/DedispersionAsm.h new file mode 100644 index 0000000000000000000000000000000000000000..d2da2718d6699c3ba97140e0c7896085a7ddcc8a --- /dev/null +++ b/RTCP/CNProc/src/DedispersionAsm.h @@ -0,0 +1,21 @@ +#ifndef LOFAR_CNPROC_DEDISPERSIONASM_H +#define LOFAR_CNPROC_DEDISPERSIONASM_H + +#if defined HAVE_BGP + +#include <Common/lofar_complex.h> + +namespace LOFAR { +namespace RTCP { + +extern "C" +{ + void _apply_chirp(fcomplex *xBuffer, fcomplex *yBuffer, const fcomplex *chirp, unsigned count); +} + +} // namespace LOFAR::RTCP +} // namespace LOFAR + +#endif + +#endif diff --git a/RTCP/CNProc/test/CMakeLists.txt b/RTCP/CNProc/test/CMakeLists.txt index d9ea03eb5f3a8889b064e2efe59ede6441fb7cf5..5e0f7013410d28f1bb402b22da63953800366ddb 100644 --- a/RTCP/CNProc/test/CMakeLists.txt +++ b/RTCP/CNProc/test/CMakeLists.txt @@ -6,6 +6,7 @@ include_directories(${PACKAGE_SOURCE_DIR}/src) lofar_add_test(tCN_Processing tCN_Processing.cc) lofar_add_test(tBeamForming tBeamForming.cc) +lofar_add_test(tDedispersion tDedispersion.cc) lofar_add_test(tPencilBeamFormer tPencilBeamFormer.cc) lofar_add_test(tStokes tStokes.cc) lofar_add_test(tInversePPF tInversePPF.cc) diff --git a/RTCP/CNProc/test/tDedispersion.cc b/RTCP/CNProc/test/tDedispersion.cc new file mode 100644 index 0000000000000000000000000000000000000000..ca4a2cbfeabd68da76ea1cf7d49a197158ffe330 --- /dev/null +++ b/RTCP/CNProc/test/tDedispersion.cc @@ -0,0 +1,121 @@ +#include <lofar_config.h> + +#include <Common/LofarLogger.h> +#include <Common/Timer.h> +#include <CN_Math.h> +#include <Dedispersion.h> + +#include <cassert> +#include <cstring> + + +#define BLOCK_SIZE 4096 +#define FFT_SIZE 4096 +#define DM 10 +#define NR_STATIONS 64 +#define NR_BEAMS 64 +#define NR_CHANNELS 16 + + +using namespace LOFAR; +using namespace LOFAR::RTCP; +using namespace LOFAR::TYPES; + + +void init(CN_Configuration &configuration) +{ + assert(BLOCK_SIZE % FFT_SIZE == 0); + + configuration.nrStations() = NR_STATIONS; + configuration.flysEye() = false; + configuration.nrPencilBeams() = NR_BEAMS; + configuration.nrChannelsPerSubband() = NR_CHANNELS; + configuration.nrSamplesPerIntegration() = BLOCK_SIZE; + configuration.dedispersionFFTsize() = FFT_SIZE; + configuration.sampleRate() = 195312.5; + configuration.dispersionMeasure() = DM; + configuration.refFreqs().push_back(50 * 195312.5); +} + + +void setTestPattern(FilteredData &filteredData) +{ + memset(&filteredData.samples[0][0][0][0], 0, filteredData.samples.num_elements() * sizeof(fcomplex)); + + for (unsigned i = 0; i < BLOCK_SIZE; i ++) + filteredData.samples[0][0][i][0] = cosisin(2 * M_PI * i * 5 / BLOCK_SIZE) /* + cosisin(2 * M_PI * i * 22 / BLOCK_SIZE) */; +} + + +void setTestPattern(BeamFormedData &beamFormedData) +{ + memset(&beamFormedData.samples[0][0][0][0], 0, beamFormedData.samples.num_elements() * sizeof(fcomplex)); + + for (unsigned i = 0; i < BLOCK_SIZE; i ++) + beamFormedData.samples[0][0][i][0] = cosisin(2 * M_PI * i * 5 / BLOCK_SIZE) /* + cosisin(2 * M_PI * i * 22 / BLOCK_SIZE) */; +} + + +void plot(const FilteredData &filteredData, float r, float g, float b) +{ + std::cout << "newcurve linetype solid linethickness 3 marktype none color " << r << ' ' << g << ' ' << b << " pts" << std::endl; + + for (unsigned i = 0; i < FFT_SIZE; i ++) + std::cout << i << ' ' << real(filteredData.samples[0][0][i][0]) << std::endl; +} + + +void plot(const BeamFormedData &beamFormedData, float r, float g, float b) +{ + std::cout << "newcurve linetype solid linethickness 3 marktype none color " << r << ' ' << g << ' ' << b << " pts" << std::endl; + + for (unsigned i = 0; i < FFT_SIZE; i ++) + std::cout << i << ' ' << real(beamFormedData.samples[0][0][i][0]) << std::endl; +} + + +int main() +{ + INIT_LOGGER_WITH_SYSINFO("tDedispersion"); + + CN_Configuration configuration; + init(configuration); + +#if 0 + BeamFormedData beamFormedData(NR_BEAMS, NR_CHANNELS, BLOCK_SIZE); + beamFormedData.allocate(); + setTestPattern(beamFormedData); + + std::cout << "newgraph xaxis size 7 yaxis size 7" << std::endl; + plot(beamFormedData, 1, 0, 0); + + std::vector<unsigned> subbands(1, 50); + DedispersionAfterBeamForming dedispersion(configuration, &beamFormedData, subbands); + + NSTimer timer("dedisperse total", true, true); + timer.start(); + dedispersion.dedisperse(&beamFormedData, 50); + timer.stop(); + + plot(beamFormedData, 0, 0, 1); +#else + FilteredData filteredData(NR_STATIONS, NR_CHANNELS, BLOCK_SIZE); + filteredData.allocate(); + setTestPattern(filteredData); + + std::cout << "newgraph xaxis size 7 yaxis size 7" << std::endl; + plot(filteredData, 1, 0, 0); + + std::vector<unsigned> subbands(1, 50); + DedispersionBeforeBeamForming dedispersion(configuration, &filteredData, subbands); + + NSTimer timer("dedisperse total", true, true); + timer.start(); + dedispersion.dedisperse(&filteredData, 50); + timer.stop(); + + plot(filteredData, 0, 0, 1); +#endif + + return 0; +} diff --git a/RTCP/Interface/include/Interface/CN_Configuration.h b/RTCP/Interface/include/Interface/CN_Configuration.h index 15b5bcf82b12d32d57fe9349b5a5edcd10f3fd3a..a6f7dcb6bdfca1fe8638635325bef3fdd124ab6f 100644 --- a/RTCP/Interface/include/Interface/CN_Configuration.h +++ b/RTCP/Interface/include/Interface/CN_Configuration.h @@ -44,6 +44,8 @@ class CN_Configuration double &startTime(); double &stopTime(); double &integrationTime(); + double &dispersionMeasure(); + unsigned &dedispersionFFTsize(); unsigned &nrStations(); unsigned nrMergedStations(); @@ -108,6 +110,8 @@ class CN_Configuration double itsStartTime; double itsStopTime; double itsIntegrationTime; + double itsDispersionMeasure; + unsigned itsDedispersionFFTsize; unsigned itsNrStations; unsigned itsNrBitsPerSample; unsigned itsNrSubbands; @@ -170,6 +174,18 @@ inline double &CN_Configuration::integrationTime() } +inline double &CN_Configuration::dispersionMeasure() +{ + return itsMarshalledData.itsDispersionMeasure; +} + + +inline unsigned &CN_Configuration::dedispersionFFTsize() +{ + return itsMarshalledData.itsDedispersionFFTsize; +} + + inline unsigned &CN_Configuration::nrStations() { return itsMarshalledData.itsNrStations; diff --git a/RTCP/Interface/include/Interface/Parset.h b/RTCP/Interface/include/Interface/Parset.h index a1d3ad48491cf62bbeb7c34be8b98028812adb23..3f3c64be755c2f17839cb8ddc01f3d97f74856b7 100644 --- a/RTCP/Interface/include/Interface/Parset.h +++ b/RTCP/Interface/include/Interface/Parset.h @@ -82,6 +82,8 @@ public: string positionType() const; vector<double> getRefPhaseCentres() const; vector<double> getPhaseCentresOf(const string &name) const; + double dispersionMeasure() const; + uint32 dedispersionFFTsize() const; uint32 CNintegrationSteps() const; uint32 IONintegrationSteps() const; uint32 stokesIntegrationSteps() const; @@ -307,6 +309,16 @@ inline double Parset::sampleDuration() const return 1.0 / sampleRate(); } +inline double Parset::dispersionMeasure() const +{ + return getDouble("OLAP.dispersionMeasure"); +} + +inline uint32 Parset::dedispersionFFTsize() const +{ + return isDefined("OLAP.CNProc.dedispersionFFTsize") ? getUint32("OLAP.CNProc.dedispersionFFTsize") : CNintegrationSteps(); +} + inline uint32 Parset::nrBitsPerSample() const { return getUint32("OLAP.nrBitsPerSample"); diff --git a/RTCP/Interface/src/CN_Configuration.cc b/RTCP/Interface/src/CN_Configuration.cc index 5086963913df9b4dec7316e0b9346c28c0745b00..88d74b910855f00896fa62a6735619c5f97c91a3 100644 --- a/RTCP/Interface/src/CN_Configuration.cc +++ b/RTCP/Interface/src/CN_Configuration.cc @@ -49,6 +49,8 @@ CN_Configuration::CN_Configuration(const Parset &parset) startTime() = parset.startTime(); stopTime() = parset.stopTime(); integrationTime() = parset.CNintegrationTime(); + dispersionMeasure() = parset.dispersionMeasure(); + dedispersionFFTsize() = parset.dedispersionFFTsize(); nrStations() = parset.nrStations(); nrBitsPerSample() = parset.nrBitsPerSample(); diff --git a/RTCP/Interface/src/Parset.cc b/RTCP/Interface/src/Parset.cc index 54faa68bc97cf500f715d3f212abdff124d890a9..3ae33fe587e8f41deef37f08adf63a3178096df8 100644 --- a/RTCP/Interface/src/Parset.cc +++ b/RTCP/Interface/src/Parset.cc @@ -122,6 +122,9 @@ void Parset::check() const { checkSubbandCount("Observation.beamList"); checkInputConsistency(); + + if (CNintegrationSteps() % dedispersionFFTsize() != 0) + THROW(InterfaceException, "OLAP.CNProc.integrationSteps (" << CNintegrationSteps() << ") must be divisible by OLAP.CNProc.dedispersionFFTsize (" << dedispersionFFTsize() << ')'); } diff --git a/RTCP/Run/src/LOFAR/Parset.py b/RTCP/Run/src/LOFAR/Parset.py index 8217a88016af302636b265a943c6f6ca76eaf33d..19e69ccdffd9c3d0fa6b1b5d699637f658b6a57a 100644 --- a/RTCP/Run/src/LOFAR/Parset.py +++ b/RTCP/Run/src/LOFAR/Parset.py @@ -138,6 +138,7 @@ class Parset(util.Parset.Parset): self.setdefault('OLAP.Storage.CoherentStokes.namemask','L${OBSID}_B${BEAM}_S${STOKES}_P${PART}_bf.raw') self.setdefault('OLAP.Storage.IncoherentStokes.namemask','L${OBSID}_SB${SUBBAND}_bf.incoherentstokes') self.setdefault('OLAP.Storage.Trigger.namemask','L${OBSID}_B${BEAM}_S${STOKES}_P${PART}_bf.trigger') + self.setdefault('OLAP.dispersionMeasure', 0); self.setdefault('OLAP.Storage.Filtered.dirmask','L${YEAR}_${OBSID}') self.setdefault('OLAP.Storage.Beamformed.dirmask','L${YEAR}_${OBSID}')