Added stokes component calculation

psrdada_cpp/effelsberg/edd/GatedSpectrometer.cuh

@@ -26,20 +26,91 @@ namespace edd {
 
 typedef unsigned long long int uint64_cu;
 static_assert(sizeof(uint64_cu) == sizeof(uint64_t), "Long long int not of 64 bit! This is problematic for CUDA!");
+
+  typedef uint64_t RawVoltageType;
+  typedef float UnpackedVoltageType;
+  typedef float2 ChannelisedVoltageType;
+
+  typedef float IntegratedPowerType;
+  //typedef int8_t IntegratedPowerType;
+
+    // Input data and intermediate processing data for one polarization
+    struct PolarizationData
+    {
+        DoubleDeviceBuffer<RawVoltageType> _raw_voltage;
+        DoubleDeviceBuffer<uint64_t> _sideChannelData;
+
+        thrust::device_vector<UnpackedVoltageType> _baseLineG0;
+        thrust::device_vector<UnpackedVoltageType> _baseLineG1;
+        thrust::device_vector<ChannelisedVoltageType> _channelised_voltage_G0;
+        thrust::device_vector<ChannelisedVoltageType> _channelised_voltage_G1;
+
+        void swap()
+        {
+            _raw_voltage.swap();
+            _sideChannelData.swap();
+        }
+    };
+
+
+    // Outptu data for one gate
+    class StokesOutput
+    {
+        public:
+      DoubleDeviceBuffer<IntegratedPowerType> I;
+      DoubleDeviceBuffer<IntegratedPowerType> Q;
+      DoubleDeviceBuffer<IntegratedPowerType> U;
+      DoubleDeviceBuffer<IntegratedPowerType> V;
+
+      DoubleDeviceBuffer<uint64_cu> _noOfBitSets;
+
+      void reset(cudaStream_t &_proc_stream)
+      {
+        thrust::fill(thrust::cuda::par.on(_proc_stream),I.a().begin(), I.a().end(), 0.);
+        thrust::fill(thrust::cuda::par.on(_proc_stream),Q.a().begin(), Q.a().end(), 0.);
+        thrust::fill(thrust::cuda::par.on(_proc_stream),U.a().begin(), U.a().end(), 0.);
+        thrust::fill(thrust::cuda::par.on(_proc_stream),V.a().begin(), V.a().end(), 0.);
+
+        thrust::fill(thrust::cuda::par.on(_proc_stream), _noOfBitSets.a().begin(), _noOfBitSets.a().end(), 0L);
+      }
+
+      void swap()
+      {
+        I.swap();
+        Q.swap();
+        U.swap();
+        V.swap();
+        _noOfBitSets.swap();
+
+      }
+
+      void resize(size_t size, size_t blocks)
+      {
+        I.resize(size * blocks);
+        Q.resize(size * blocks);
+        U.resize(size * blocks);
+        V.resize(size * blocks);
+
+        _noOfBitSets.resize(blocks);
+      }
+    };
+
+
+
+
+
+
+
 /**
  @class GatedSpectrometer
  @brief Split data into two streams and create integrated spectra depending on
  bit set in side channel data.
 
  */
-template <class HandlerType, typename IntegratedPowerType> class GatedSpectrometer {
+template <class HandlerType> class GatedSpectrometer {
 public:
-  typedef uint64_t RawVoltageType;
-  typedef float UnpackedVoltageType;
-  typedef float2 ChannelisedVoltageType;
 
-//  typedef float IntegratedPowerType;
-  //typedef int8_t IntegratedPowerType;
+
 
 public:
   /**
@@ -90,11 +161,10 @@ public:
   bool operator()(RawBytes &block);
 
 private:
-  void process(thrust::device_vector<RawVoltageType> const &digitiser_raw,
-               thrust::device_vector<uint64_t> const &sideChannelData,
-               thrust::device_vector<IntegratedPowerType> &detected,
-               thrust::device_vector<uint64_cu> &noOfBitSetsIn_G0,
-               thrust::device_vector<uint64_cu> &noOfBitSetsIn_G1);
+  // gate the data and fft data per gate
+  void gated_fft(PolarizationData &data,
+  thrust::device_vector<uint64_cu> &_noOfBitSetsIn_G0,
+  thrust::device_vector<uint64_cu> &_noOfBitSetsIn_G1);
 
 private:
   DadaBufferLayout _dadaBufferLayout;
@@ -115,26 +185,20 @@ private:
   double _processing_efficiency;
 
   std::unique_ptr<Unpacker> _unpacker;
-  std::unique_ptr<DetectorAccumulator<IntegratedPowerType> > _detector;
 
-  // Input data
-  DoubleDeviceBuffer<RawVoltageType> _raw_voltage_db;
-  DoubleDeviceBuffer<uint64_t> _sideChannelData_db;
+  // Input data and per pol intermediate data
+  PolarizationData polarization0, polarization1;
 
   // Output data
-  DoubleDeviceBuffer<IntegratedPowerType> _power_db;
-
+  StokesOutput stokes_G0, stokes_G1;
 
-  DoubleDeviceBuffer<uint64_cu> _noOfBitSetsIn_G0;
-  DoubleDeviceBuffer<uint64_cu> _noOfBitSetsIn_G1;
   DoublePinnedHostBuffer<char> _host_power_db;
 
-  // Intermediate process steps
+  // Temporary processing block
+  // ToDo: Use inplace FFT to avoid temporary coltage array
   thrust::device_vector<UnpackedVoltageType> _unpacked_voltage_G0;
   thrust::device_vector<UnpackedVoltageType> _unpacked_voltage_G1;
-  thrust::device_vector<ChannelisedVoltageType> _channelised_voltage;
-  thrust::device_vector<UnpackedVoltageType> _baseLineNG0;
-  thrust::device_vector<UnpackedVoltageType> _baseLineNG1;
+
 
   cudaStream_t _h2d_stream;
   cudaStream_t _proc_stream;
@@ -142,11 +206,10 @@ private:
 };
 
 
-
 /**
    * @brief      Splits the input data depending on a bit set into two arrays.
    *
-   * @detail     The resulting gaps are filled with zeros in the other stream.
+   * @detail     The resulting gaps are filled with a given baseline value in the other stream.
    *
    * @param      GO Input data. Data is set to the baseline value if corresponding
    *             sideChannelData bit at bitpos os set.
@@ -177,6 +240,57 @@ __global__ void gating(float *G0, float *G1, const int64_t *sideChannelData,
                        uint64_cu* stats_G0,
                        uint64_cu* stats_G1);
 
+/**
+ * @brief calculate stokes IQUV from two complex valuies for each polarization
+ */
+//__host__ __device__ void stokes_IQUV(const float2 &p1, const float2 &p2, float &I, float &Q, float &U, float &V);
+__host__ __device__ void stokes_IQUV(const float2 &p1, const float2 &p2, float &I, float &Q, float &U, float &V)
+{
+    I = fabs(p1.x*p1.x + p1.y * p1.y) + fabs(p2.x*p2.x + p2.y * p2.y);
+    Q = fabs(p1.x*p1.x + p1.y * p1.y) - fabs(p2.x*p2.x + p2.y * p2.y);
+    U = 2 * (p1.x*p2.x + p1.y * p2.y);
+    V = -2 * (p1.y*p2.x - p1.x * p2.y);
+}
+
+
+
+
+/**
+ * @brief calculate stokes IQUV spectra pol1, pol2 are arrays of naccumulate
+ * complex spectra for individual polarizations
+ */
+__global__ void stokes_accumulate(float2 const __restrict__ *pol1,
+        float2 const __restrict__ *pol2, float *I, float* Q, float *U, float*V,
+        int nchans, int naccumulate)
+{
+
+  for (size_t i = blockIdx.x * blockDim.x + threadIdx.x; (i < nchans);
+       i += blockDim.x * gridDim.x)
+  {
+      float rI = 0;
+      float rQ = 0;
+      float rU = 0;
+      float rV = 0;
+
+      for (int k=0; k < naccumulate; k++)
+      {
+        const float2 p1 = pol1[i + k * nchans];
+        const float2 p2 = pol2[i + k * nchans];
+
+        rI += fabs(p1.x * p1.x + p1.y * p1.y) + fabs(p2.x * p2.x + p2.y * p2.y);
+        rQ += fabs(p1.x * p1.x + p1.y * p1.y) - fabs(p2.x * p2.x + p2.y * p2.y);
+        rU += 2.f * (p1.x * p2.x + p1.y * p2.y);
+        rV += -2.f * (p1.y * p2.x - p1.x * p2.y);
+      }
+      I[i] += rI;
+      Q[i] += rQ;
+      U[i] += rU;
+      V[i] += rV;
+  }
+
+}
+
+
 
 
 

psrdada_cpp/effelsberg/edd/src/GatedSpectrometer_cli.cu

@@ -33,7 +33,7 @@ void launchSpectrometer(const effelsberg::edd::DadaBufferLayout &dadaBufferLayou
     if (output_type == "file")
     {
       SimpleFileWriter sink(filename);
-      effelsberg::edd::GatedSpectrometer<decltype(sink), T> spectrometer(dadaBufferLayout,
+      effelsberg::edd::GatedSpectrometer<decltype(sink)> spectrometer(dadaBufferLayout,
           selectedSideChannel, selectedBit,
           fft_length, naccumulate, nbits, input_level,
           output_level, sink);
@@ -45,7 +45,7 @@ void launchSpectrometer(const effelsberg::edd::DadaBufferLayout &dadaBufferLayou
     else if (output_type == "dada")
     {
       DadaOutputStream sink(string_to_key(filename), log);
-      effelsberg::edd::GatedSpectrometer<decltype(sink), T> spectrometer(dadaBufferLayout,
+      effelsberg::edd::GatedSpectrometer<decltype(sink)> spectrometer(dadaBufferLayout,
           selectedSideChannel, selectedBit,
           fft_length, naccumulate, nbits, input_level,
           output_level, sink);
@@ -185,13 +185,8 @@ int main(int argc, char **argv) {
 
     effelsberg::edd::DadaBufferLayout bufferLayout(input_key, speadHeapSize, nSideChannels);
 
-    if (output_bit_depth == 8)
-    {
-      launchSpectrometer<int8_t>(bufferLayout, output_type, filename,
-          selectedSideChannel, selectedBit,
-       fft_length, naccumulate, nbits, input_level, output_level);
-    }
-    else if (output_bit_depth == 32)
+    // ToDo: Supprot only single output depth
+    if (output_bit_depth == 32)
     {
       launchSpectrometer<float>(bufferLayout, output_type, filename,
           selectedSideChannel, selectedBit,

psrdada_cpp/effelsberg/edd/test/src/GatedSpectrometerTest.cu

@@ -7,7 +7,6 @@
 #include "thrust/device_vector.h"
 #include "thrust/extrema.h"
 
-namespace {
 
 TEST(GatedSpectrometer, BitManipulationMacros) {
   for (int i = 0; i < 64; i++) {
@@ -23,31 +22,123 @@ TEST(GatedSpectrometer, BitManipulationMacros) {
   }
 }
 
-//
-//TEST(GatedSpectrometer, ParameterSanity) {
-//  ::testing::FLAGS_gtest_death_test_style = "threadsafe";
-//  psrdada_cpp::NullSink sink;
-//
-//  // 8 or 12 bit sampling
-//  EXPECT_DEATH(psrdada_cpp::effelsberg::edd::GatedSpectrometer<decltype(sink),int8_t > (0, 0, 0, 0, 4096, 0, 0, 0, 0, 0, sink),
-//               "_nbits == 8");
-//  // naccumulate > 0
-//  EXPECT_DEATH(psrdada_cpp::effelsberg::edd::GatedSpectrometer<decltype(sink),int8_t > (0, 0, 0, 0, 4096, 0, 0, 8, 0, 0, sink),
-//               "_naccumulate");
-//
-//  // selected side channel
-//  EXPECT_DEATH(psrdada_cpp::effelsberg::edd::GatedSpectrometer<decltype(sink),int8_t > (0, 1, 2, 0, 4096, 0, 1, 8, 0, 0, sink),
-//               "nSideChannels");
-//
-//  // selected bit
-//  EXPECT_DEATH(psrdada_cpp::effelsberg::edd::GatedSpectrometer<decltype(sink),int8_t > (0, 2, 1, 65, 4096, 0, 1, 8, 0, 0, sink),
-//               "selectedBit");
-//
-//  // valid construction
-//  psrdada_cpp::effelsberg::edd::GatedSpectrometer<decltype(sink), int8_t> a(
-//      4096 * 4096, 2, 1, 63, 4096, 1024, 1, 8, 100., 100., sink);
-//}
-} // namespace
+
+TEST(GatedSpectrometer, stokes_IQUV)
+{
+    float I,Q,U,V;
+    // No field
+    psrdada_cpp::effelsberg::edd::stokes_IQUV((float2){0.0f,0.0f}, (float2){0.0f,0.0f}, I, Q, U, V);
+    EXPECT_FLOAT_EQ(I, 0);
+    EXPECT_FLOAT_EQ(Q, 0);
+    EXPECT_FLOAT_EQ(U, 0);
+    EXPECT_FLOAT_EQ(V, 0);
+
+    // For p1 = Ex, p2 = Ey
+    // horizontal polarized
+    psrdada_cpp::effelsberg::edd::stokes_IQUV((float2){1.0f,0.0f}, (float2){0.0f,0.0f}, I, Q, U, V);
+    EXPECT_FLOAT_EQ(I, 1);
+    EXPECT_FLOAT_EQ(Q, 1);
+    EXPECT_FLOAT_EQ(U, 0);
+    EXPECT_FLOAT_EQ(V, 0);
+
+
+    // vertical polarized
+    psrdada_cpp::effelsberg::edd::stokes_IQUV((float2){0.0f,0.0f}, (float2){1.0f,0.0f}, I, Q, U, V);
+    EXPECT_FLOAT_EQ(I, 1);
+    EXPECT_FLOAT_EQ(Q, -1);
+    EXPECT_FLOAT_EQ(U, 0);
+    EXPECT_FLOAT_EQ(V, 0);
+
+    //linear +45 deg.
+    psrdada_cpp::effelsberg::edd::stokes_IQUV((float2){1.0f/std::sqrt(2),0.0f}, (float2){1.0f/std::sqrt(2),0.0f}, I, Q, U, V);
+    EXPECT_FLOAT_EQ(I, 1);
+    EXPECT_FLOAT_EQ(Q, 0);
+    EXPECT_FLOAT_EQ(U, 1);
+    EXPECT_FLOAT_EQ(V, 0);
+
+    //linear -45 deg.
+    psrdada_cpp::effelsberg::edd::stokes_IQUV((float2){-1.0f/std::sqrt(2),0.0f}, (float2){1.0f/std::sqrt(2),0.0f}, I, Q, U, V);
+    EXPECT_FLOAT_EQ(I, 1);
+    EXPECT_FLOAT_EQ(Q, 0);
+    EXPECT_FLOAT_EQ(U, -1);
+    EXPECT_FLOAT_EQ(V, 0);
+
+
+
+    //left circular
+    psrdada_cpp::effelsberg::edd::stokes_IQUV((float2){.0f,1.0f/std::sqrt(2)}, (float2){1.0f/std::sqrt(2),.0f}, I, Q, U, V);
+    EXPECT_FLOAT_EQ(I, 1);
+    EXPECT_FLOAT_EQ(Q, 0);
+    EXPECT_FLOAT_EQ(U, 0);
+    EXPECT_FLOAT_EQ(V, -1);
+
+    // right circular
+    psrdada_cpp::effelsberg::edd::stokes_IQUV((float2){.0f,-1.0f/std::sqrt(2)}, (float2){1.0f/std::sqrt(2),.0f}, I, Q, U, V);
+    EXPECT_FLOAT_EQ(I, 1);
+    EXPECT_FLOAT_EQ(Q, 0);
+    EXPECT_FLOAT_EQ(U, 0);
+    EXPECT_FLOAT_EQ(V, 1);
+
+}
+
+
+TEST(GatedSpectrometer, stokes_accumulate)
+{
+    size_t nchans = 8 * 1024 * 1024 + 1;
+    size_t naccumulate = 5;
+
+    thrust::device_vector<float2> P0(nchans * naccumulate);
+    thrust::device_vector<float2> P1(nchans * naccumulate);
+    thrust::fill(P0.begin(), P0.end(), (float2){0, 0});
+    thrust::fill(P1.begin(), P1.end(), (float2){0, 0});
+    thrust::device_vector<float> I(nchans);
+    thrust::device_vector<float> Q(nchans);
+    thrust::device_vector<float> U(nchans);
+    thrust::device_vector<float> V(nchans);
+    thrust::fill(I.begin(), I.end(), 0);
+    thrust::fill(Q.begin(), Q.end(), 0);
+    thrust::fill(U.begin(), U.end(), 0);
+    thrust::fill(V.begin(), V.end(), 0);
+
+    // This channel should be left circular polarized
+    size_t idx0 = 23;
+    for (int k = 0; k< naccumulate; k++)
+    {
+        size_t idx = idx0 + k * nchans;
+        P0[idx] = (float2){0.0f, 1.0f/std::sqrt(2)};
+        P1[idx] = (float2){1.0f/std::sqrt(2),0.0f};
+    }
+
+    psrdada_cpp::effelsberg::edd::stokes_accumulate<<<1024, 1024>>>(
+          thrust::raw_pointer_cast(P0.data()),
+          thrust::raw_pointer_cast(P1.data()),
+          thrust::raw_pointer_cast(I.data()),
+          thrust::raw_pointer_cast(Q.data()),
+          thrust::raw_pointer_cast(U.data()),
+          thrust::raw_pointer_cast(V.data()),
+          nchans,
+          naccumulate
+            );
+
+    thrust::pair<thrust::device_vector<float>::iterator, thrust::device_vector<float>::iterator> minmax;
+
+    minmax = thrust::minmax_element(I.begin(), I.end());
+    EXPECT_FLOAT_EQ(*minmax.first, 0);
+    EXPECT_FLOAT_EQ(*minmax.second, naccumulate);
+
+    minmax = thrust::minmax_element(Q.begin(), Q.end());
+    EXPECT_FLOAT_EQ(*minmax.first, 0);
+    EXPECT_FLOAT_EQ(*minmax.second, 0);
+
+    minmax = thrust::minmax_element(U.begin(), U.end());
+    EXPECT_FLOAT_EQ(*minmax.first, 0);
+    EXPECT_FLOAT_EQ(*minmax.second, 0);
+
+    minmax = thrust::minmax_element(V.begin(), V.end());
+    EXPECT_FLOAT_EQ(*minmax.first, -1. * naccumulate);
+    EXPECT_FLOAT_EQ(*minmax.second, 0);
+
+}
 
 
 TEST(GatedSpectrometer, GatingKernel)