From 0ac466e55ded15d3b00ad4b776a0c15c11042eba Mon Sep 17 00:00:00 2001
From: Ewan Barr <ewanbarr@krpc74.mpifr-bonn.mpg.de>
Date: Tue, 11 Dec 2018 13:15:42 +0100
Subject: [PATCH] changed offset and scaling calculation

---
 psrdada_cpp/effelsberg/edd/FftSpectrometer.cuh      |  3 +--
 .../effelsberg/edd/detail/FftSpectrometer.cu        | 13 ++++++++-----
 .../effelsberg/edd/src/fft_spectrometer_cli.cu      | 11 +++++------
 .../edd/test/src/FftSpectrometerTester.cu           |  2 +-
 4 files changed, 15 insertions(+), 14 deletions(-)

diff --git a/psrdada_cpp/effelsberg/edd/FftSpectrometer.cuh b/psrdada_cpp/effelsberg/edd/FftSpectrometer.cuh
index 75202028..ca3a6862 100644
--- a/psrdada_cpp/effelsberg/edd/FftSpectrometer.cuh
+++ b/psrdada_cpp/effelsberg/edd/FftSpectrometer.cuh
@@ -63,8 +63,7 @@ private:
     std::size_t _fft_length;
     std::size_t _naccumulate;
     std::size_t _nbits;
-    float _scaling;
-    float _offset;
+    float _input_level;
     HandlerType& _handler;
     cufftHandle _fft_plan;
     int _nchans;
diff --git a/psrdada_cpp/effelsberg/edd/detail/FftSpectrometer.cu b/psrdada_cpp/effelsberg/edd/detail/FftSpectrometer.cu
index edce0e13..33d464a9 100644
--- a/psrdada_cpp/effelsberg/edd/detail/FftSpectrometer.cu
+++ b/psrdada_cpp/effelsberg/edd/detail/FftSpectrometer.cu
@@ -14,15 +14,13 @@ FftSpectrometer<HandlerType>::FftSpectrometer(
     std::size_t fft_length,
     std::size_t naccumulate,
     std::size_t nbits,
-    float scaling,
-    float offset,
+    float input_level,
     HandlerType& handler)
     : _buffer_bytes(buffer_bytes)
     , _fft_length(fft_length)
     , _naccumulate(naccumulate)
     , _nbits(nbits)
-    , _scaling(scaling)
-    , _offset(offset)
+    , _input_level(input_level)
     , _handler(handler)
     , _fft_plan(0)
     , _call_count(0)
@@ -38,6 +36,11 @@ FftSpectrometer<HandlerType>::FftSpectrometer(
     std::size_t n64bit_words = buffer_bytes / sizeof(uint64_t);
     _nchans = _fft_length / 2 + 1;
     int batch = nsamps_per_buffer/_fft_length;
+    BOOST_LOG_TRIVIAL(debug) << "Calculating scales and offsets";
+    float dof = 2 * _naccumulate;
+    float scale = (_input_level * np.sqrt(_nchans))**2;
+    float offset = scale * dof;
+    float scaling = scale * std::sqrt(2 * dof);
     BOOST_LOG_TRIVIAL(debug) << "Generating FFT plan";
     int n[] = {static_cast<int>(_fft_length)};
     CUFFT_ERROR_CHECK(cufftPlanMany(&_fft_plan, 1, n, NULL, 1, _fft_length,
@@ -55,7 +58,7 @@ FftSpectrometer<HandlerType>::FftSpectrometer(
     CUFFT_ERROR_CHECK(cufftSetStream(_fft_plan, _proc_stream));
     _unpacker.reset(new Unpacker(_proc_stream));
     _detector.reset(new DetectorAccumulator(_nchans, _naccumulate,
-        _scaling, _offset, _proc_stream));
+        scaling, offset, _proc_stream));
 }
 
 template <class HandlerType>
diff --git a/psrdada_cpp/effelsberg/edd/src/fft_spectrometer_cli.cu b/psrdada_cpp/effelsberg/edd/src/fft_spectrometer_cli.cu
index 537108e8..f929813a 100644
--- a/psrdada_cpp/effelsberg/edd/src/fft_spectrometer_cli.cu
+++ b/psrdada_cpp/effelsberg/edd/src/fft_spectrometer_cli.cu
@@ -30,8 +30,7 @@ int main(int argc, char** argv)
         int nsamps_per_block;
         int naccumulate;
         int nbits;
-        float scaling;
-        float offset;
+        float input_level;
         std::time_t now = std::time(NULL);
         std::tm * ptm = std::localtime(&now);
         char buffer[32];
@@ -62,8 +61,8 @@ int main(int argc, char** argv)
         ("nbits,b", po::value<int>(&nbits)->required(),
             "The number of bits per sample in the packetiser output (8 or 12)")
 
-        ("scaling", po::value<float>(&scaling)->required(),
-            "The power scaling for data produced by the spectrometer (used for conversion back to 8-bit)")
+        ("input_level", po::value<float>(&input_level)->required(),
+            "The input power level (standard deviation, used for 8-bit conversion)")
 
         ("offset", po::value<float>(&offset)->required(),
             "The power offset for data produced by the spectrometer (used for conversion back to 8-bit)")
@@ -86,7 +85,7 @@ int main(int argc, char** argv)
             po::store(po::parse_command_line(argc, argv, desc), vm);
             if ( vm.count("help")  )
             {
-                std::cout << "EDDFFT -- Read EDD data from a DADA buffer and perform a simple FFT spectrometer"
+                std::cout << "FftSpectrometer -- Read EDD data from a DADA buffer and perform a simple FFT spectrometer"
                 << std::endl << desc << std::endl;
                 return SUCCESS;
             }
@@ -106,7 +105,7 @@ int main(int argc, char** argv)
         std::size_t buffer_bytes = client.data_buffer_size()
         SimpleFileWriter sink(filename);
         //NullSink sink;
-        effelsberg::edd::FftSpectrometer<decltype(sink)> spectrometer(buffer_bytes, fft_length, naccumulate, nbits, scaling, offset, sink);
+        effelsberg::edd::FftSpectrometer<decltype(sink)> spectrometer(buffer_bytes, fft_length, naccumulate, nbits, input_level, sink);
         DadaInputStream<decltype(spectrometer)> istream(input_key, log, spectrometer);
         istream.start();
         /**
diff --git a/psrdada_cpp/effelsberg/edd/test/src/FftSpectrometerTester.cu b/psrdada_cpp/effelsberg/edd/test/src/FftSpectrometerTester.cu
index 15929fae..efb8b709 100644
--- a/psrdada_cpp/effelsberg/edd/test/src/FftSpectrometerTester.cu
+++ b/psrdada_cpp/effelsberg/edd/test/src/FftSpectrometerTester.cu
@@ -42,7 +42,7 @@ void FftSpectrometerTester::performance_test(
     std::vector<char> header_block(4096);
     RawBytes header_raw_bytes(header_block.data(), 4096, 4096);
     NullSink null_sink;
-    FftSpectrometer<NullSink> spectrometer(input_block_bytes, fft_length, tscrunch, nbits, 1.0f, 0.0f, null_sink);
+    FftSpectrometer<NullSink> spectrometer(input_block_bytes, fft_length, tscrunch, nbits, 16.0f, null_sink);
     spectrometer.init(header_raw_bytes);
     for (int ii = 0; ii < 100; ++ii)
     {
-- 
GitLab