Commit 418e19c7 authored by Tobias Winchen's avatar Tobias Winchen
Browse files

Count saturated samples in dual spec output

parent 72008d9d
......@@ -56,6 +56,7 @@ struct PolarizationData
DoubleDeviceBuffer<RawVoltageType> _raw_voltage;
/// Side channel data
DoubleDeviceBuffer<uint64_t> _sideChannelData;
DoubleHostBuffer<uint64_t> _sideChannelData_h;
/// Baseline in gate 0 state
thrust::device_vector<UnpackedVoltageType> _baseLineG0;
......@@ -164,6 +165,10 @@ struct PowerSpectrumOutput
/// Number of samples integrated in this output block
DoubleDeviceBuffer<uint64_cu> _noOfBitSets;
/// Number of samples overflowed
DoubleHostBuffer<uint64_cu> _noOfOverflowed;
/// Swap output buffers and reset the buffer in given stream for new integration
void swap(cudaStream_t &_proc_stream);
};
......@@ -237,6 +242,9 @@ struct FullStokesOutput
/// Number of samples integrated in this output block
DoubleDeviceBuffer<uint64_cu> _noOfBitSets;
/// Number of samples overflowed
DoubleHostBuffer<uint64_cu> _noOfOverflowed;
/// Swap output buffers
void swap(cudaStream_t &_proc_stream);
};
......@@ -316,23 +324,27 @@ public:
*/
bool operator()(RawBytes &block);
private:
// Processing strategy for single pol mode
// make the relevant processing methods proteceed only for testing
protected:
/// Processing strategy for single pol mode
void process(SinglePolarizationInput *inputDataStream, GatedPowerSpectrumOutput *outputDataStream);
// Processing strategy for dual pol power mode
/// Processing strategy for dual pol power mode
//void process(DualPolarizationInput*inputDataStream, GatedPowerSpectrumOutput *outputDataStream);
// Processing strategy for full stokes mode
/// Processing strategy for full stokes mode
void process(DualPolarizationInput *inputDataStream, GatedFullStokesOutput *outputDataStream);
OutputType* outputDataStream;
InputType* inputDataStream;
private:
// gate the data from the input stream and fft data per gate. Number of bit
// sets in every gate is stored in the output datasets
void gated_fft(PolarizationData &data,
thrust::device_vector<uint64_cu> &_noOfBitSetsIn_G0,
thrust::device_vector<uint64_cu> &_noOfBitSetsIn_G1);
private:
DadaBufferLayout _dadaBufferLayout;
std::size_t _fft_length;
std::size_t _naccumulate;
......@@ -349,9 +361,6 @@ private:
std::unique_ptr<Unpacker> _unpacker;
OutputType* outputDataStream;
InputType* inputDataStream;
// Temporary processing block
// ToDo: Use inplace FFT to avoid temporary voltage array
thrust::device_vector<UnpackedVoltageType> _unpacked_voltage_G0;
......
......@@ -364,8 +364,6 @@ void GatedSpectrometer<HandlerType, InputType, OutputType>::process(SinglePolari
{
gated_fft(*inputDataStream, outputDataStream->G0._noOfBitSets.a(), outputDataStream->G1._noOfBitSets.a());
kernels::detect_and_accumulate<IntegratedPowerType> <<<1024, 1024, 0, _proc_stream>>>(
thrust::raw_pointer_cast(inputDataStream->_channelised_voltage_G0.data()),
thrust::raw_pointer_cast(outputDataStream->G0.data.a().data()),
......@@ -382,6 +380,34 @@ void GatedSpectrometer<HandlerType, InputType, OutputType>::process(SinglePolari
_naccumulate / _nBlocks,
1, 0., 1, 0);
// count saturated samples
for(int output_block_number = 0; output_block_number <_nBlocks; output_block_number++)
{
outputDataStream->G0._noOfOverflowed.a()[output_block_number] = 0;
outputDataStream->G1._noOfOverflowed.a()[output_block_number] = 0;
const int heaps_per_output_spectra = inputDataStream->_sideChannelData_h.size() / _naccumulate / _nBlocks;
for (int j = output_block_number * heaps_per_output_spectra ; j < (output_block_number+1) * heaps_per_output_spectra * _dadaBufferLayout.getNSideChannels(); j+=_dadaBufferLayout.getNSideChannels())
{
if (TEST_BIT(inputDataStream->_sideChannelData_h.a().data()[j], 42))
{ // heap was lost
continue;
}
uint16_t n_saturated = (inputDataStream->_sideChannelData_h.a().data()[j] >> 32);
if (TEST_BIT(inputDataStream->_sideChannelData_h.a().data()[j], _selectedBit))
{
outputDataStream->G1._noOfOverflowed.a()[output_block_number] += n_saturated;
}
else
{
outputDataStream->G0._noOfOverflowed.a()[output_block_number] += n_saturated;
}
}
BOOST_LOG_TRIVIAL(debug) << "Number of saturated samples G0: " << outputDataStream->G0._noOfOverflowed.a()[output_block_number]
<< ", G1: " << outputDataStream->G1._noOfOverflowed.a()[output_block_number];
}
}
......@@ -418,6 +444,43 @@ void GatedSpectrometer<HandlerType, InputType, OutputType>::process(DualPolariza
thrust::raw_pointer_cast(outputDataStream->G1.V.a().data() + output_offset),
_nchans, _naccumulate / _nBlocks
);
// count saturated samples
outputDataStream->G0._noOfOverflowed.a()[output_block_number] = 0;
outputDataStream->G1._noOfOverflowed.a()[output_block_number] = 0;
size_t lostHeaps = 0;
const int heaps_per_output_spectra = inputDataStream->polarization0._sideChannelData_h.size() / _naccumulate / _nBlocks;
for (int j = output_block_number * heaps_per_output_spectra ; j < (output_block_number+1) * heaps_per_output_spectra * _dadaBufferLayout.getNSideChannels(); j+=_dadaBufferLayout.getNSideChannels())
{
if (TEST_BIT(inputDataStream->polarization0._sideChannelData_h.a().data()[j], 42) || TEST_BIT(inputDataStream->polarization1._sideChannelData_h.a().data()[j], 42))
{
lostHeaps++;
// thiss heap was lost, do not count
continue;
}
uint16_t n_saturated = (inputDataStream->polarization0._sideChannelData_h.a().data()[j] >> 32) + (inputDataStream->polarization1._sideChannelData_h.a().data()[j] >> 32);
if (TEST_BIT(inputDataStream->polarization0._sideChannelData_h.a().data()[j], _selectedBit))
{
outputDataStream->G1._noOfOverflowed.a()[output_block_number] += n_saturated;
}
else
{
outputDataStream->G0._noOfOverflowed.a()[output_block_number] += n_saturated;
}
}
BOOST_LOG_TRIVIAL(debug) << "Number of saturated samples G0: " << outputDataStream->G0._noOfOverflowed.a()[output_block_number]
<< ", G1: " << outputDataStream->G1._noOfOverflowed.a()[output_block_number];
BOOST_LOG_TRIVIAL(debug) << "Lost heaps: " << lostHeaps;
}
// cudaDeviceSynchronize();
......
......@@ -28,7 +28,9 @@ __global__ void gating(float* __restrict__ G0,
const float* __restrict__ _baseLineG1,
float* __restrict__ baseLineNG0,
float* __restrict__ baseLineNG1,
uint64_cu* stats_G0, uint64_cu* stats_G1) {
uint64_cu* stats_G0,
uint64_cu* stats_G1)
{
// statistics values for samopels to G0, G1
uint32_t _G0stats = 0;
uint32_t _G1stats = 0;
......@@ -91,6 +93,9 @@ __global__ void gating(float* __restrict__ G0,
// Updates the baselines of the gates for the polarization set for the next
// block
// only few output blocks per input block thus execution on only one thread.
......@@ -188,6 +193,7 @@ void PolarizationData::resize(size_t rawVolttageBufferBytes, size_t nsidechannel
_channelised_voltage_G0.resize(channelized_samples);
_channelised_voltage_G1.resize(channelized_samples);
_sideChannelData.resize(nsidechannelitems);
_sideChannelData_h.resize(nsidechannelitems);
BOOST_LOG_TRIVIAL(debug) << " Channelised voltages size: " << _channelised_voltage_G0.size();
}
......@@ -196,8 +202,8 @@ SinglePolarizationInput::SinglePolarizationInput(size_t fft_length, size_t nbits
&dadaBufferLayout) : PolarizationData(nbits), _fft_length(fft_length), _dadaBufferLayout(dadaBufferLayout)
{
size_t nsamps_per_buffer = _dadaBufferLayout.sizeOfData() * 8 / nbits;
size_t _batch = nsamps_per_buffer / _fft_length;
size_t nsamps_per_buffer = _dadaBufferLayout.sizeOfData() * 8 / nbits;
size_t _batch = nsamps_per_buffer / _fft_length;
resize(_dadaBufferLayout.sizeOfData(), _dadaBufferLayout.getNSideChannels() * _dadaBufferLayout.getNHeaps(), (_fft_length / 2 + 1) * _batch);
};
......@@ -229,6 +235,11 @@ void SinglePolarizationInput::getFromBlock(RawBytes &block, cudaStream_t &_h2d_s
static_cast<void *>(_sideChannelData.a_ptr()),
static_cast<void *>(block.ptr() + _dadaBufferLayout.sizeOfData() + _dadaBufferLayout.sizeOfGap()),
_dadaBufferLayout.sizeOfSideChannelData(), cudaMemcpyHostToDevice, _h2d_stream));
std::memcpy(static_cast<void *>(_sideChannelData_h.a_ptr()),
static_cast<void *>(block.ptr() + _dadaBufferLayout.sizeOfData() + _dadaBufferLayout.sizeOfGap()),
_dadaBufferLayout.sizeOfSideChannelData());
BOOST_LOG_TRIVIAL(debug) << "First side channel item: 0x" << std::setw(16)
<< std::setfill('0') << std::hex <<
(reinterpret_cast<uint64_t*>(block.ptr() + _dadaBufferLayout.sizeOfData()
......@@ -298,6 +309,7 @@ CUDA_ERROR_CHECK(cudaMemcpy2DAsync(
_dadaBufferLayout.sizeOfSideChannelData() / 2 / sizeof(uint64_t),
cudaMemcpyHostToDevice, _h2d_stream));
CUDA_ERROR_CHECK(cudaMemcpy2DAsync(
static_cast<void *>(polarization1._sideChannelData.a_ptr()),
sizeof(uint64_t),
......@@ -306,6 +318,14 @@ CUDA_ERROR_CHECK(cudaMemcpy2DAsync(
sizeof(uint64_t),
_dadaBufferLayout.sizeOfSideChannelData() / 2 / sizeof(uint64_t), cudaMemcpyHostToDevice, _h2d_stream));
uint64_t * sci = reinterpret_cast<uint64_t *>(block.ptr() + _dadaBufferLayout.sizeOfData() + _dadaBufferLayout.sizeOfGap());
for(int i =0; i < _dadaBufferLayout.sizeOfSideChannelData() / sizeof(uint64_t); i+=2)
{
polarization0._sideChannelData_h.a()[i / 2] = sci[i];
polarization1._sideChannelData_h.a()[i / 2] = sci[i + 1];
}
BOOST_LOG_TRIVIAL(debug) << "First side channel item: 0x" << std::setw(16)
<< std::setfill('0') << std::hex <<
(reinterpret_cast<uint64_t*>(block.ptr() + _dadaBufferLayout.sizeOfData()
......@@ -319,6 +339,7 @@ PowerSpectrumOutput::PowerSpectrumOutput(size_t size, size_t blocks)
BOOST_LOG_TRIVIAL(debug) << "Setting size of power spectrum output size = " << size << ", blocks = " << blocks;
data.resize(size * blocks);
_noOfBitSets.resize(blocks);
_noOfOverflowed.resize(blocks);
}
......@@ -326,6 +347,7 @@ void PowerSpectrumOutput::swap(cudaStream_t &_proc_stream)
{
data.swap();
_noOfBitSets.swap();
_noOfOverflowed.swap();
thrust::fill(thrust::cuda::par.on(_proc_stream), data.a().begin(), data.a().end(), 0.);
thrust::fill(thrust::cuda::par.on(_proc_stream), _noOfBitSets.a().begin(), _noOfBitSets.a().end(), 0L);
}
......@@ -337,7 +359,7 @@ G1(nchans, blocks)
{
// on the host both power are stored in the same data buffer together with
// the number of bit sets
_host_power.resize( 2 * ( _nchans * sizeof(IntegratedPowerType) + sizeof(size_t) ) * G0._noOfBitSets.size());
_host_power.resize( 2 * ( _nchans * sizeof(IntegratedPowerType) + 2 * sizeof(size_t) ) * G0._noOfBitSets.size());
}
......@@ -360,7 +382,7 @@ void GatedPowerSpectrumOutput::data2Host(cudaStream_t &_d2h_stream)
// size of individual spectrum + meta
size_t memslicesize = (_nchans * sizeof(IntegratedPowerType));
// number of spectra per output
size_t memOffset = 2 * i * (memslicesize + sizeof(size_t));
size_t memOffset = 2 * i * (memslicesize + 2*sizeof(size_t));
// copy 2x channel data
CUDA_ERROR_CHECK(
......@@ -382,11 +404,20 @@ void GatedPowerSpectrumOutput::data2Host(cudaStream_t &_d2h_stream)
1 * sizeof(size_t),
cudaMemcpyDeviceToHost, _d2h_stream));
std::memcpy(static_cast<void *>(_host_power.a_ptr() + memOffset + 2 * _nchans * sizeof(IntegratedPowerType) + sizeof(size_t)),
static_cast<void *>(G0._noOfOverflowed.b_ptr() + i),
1 * sizeof(size_t));
CUDA_ERROR_CHECK(
cudaMemcpyAsync( static_cast<void *>(_host_power.a_ptr() + memOffset + 2 * _nchans * sizeof(IntegratedPowerType) + sizeof(size_t)),
cudaMemcpyAsync( static_cast<void *>(_host_power.a_ptr() + memOffset + 2 * _nchans * sizeof(IntegratedPowerType) + 2 * sizeof(size_t)),
static_cast<void *>(G1._noOfBitSets.b_ptr() + i ),
1 * sizeof(size_t),
cudaMemcpyDeviceToHost, _d2h_stream));
std::memcpy(static_cast<void *>(_host_power.a_ptr() + memOffset + 2 * _nchans * sizeof(IntegratedPowerType) + 3*sizeof(size_t)),
static_cast<void *>(G1._noOfOverflowed.b_ptr() + i),
1 * sizeof(size_t));
}
}
......@@ -398,6 +429,7 @@ void FullStokesOutput::swap(cudaStream_t &_proc_stream)
U.swap();
V.swap();
_noOfBitSets.swap();
_noOfOverflowed.swap();
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.);
......@@ -413,6 +445,7 @@ FullStokesOutput::FullStokesOutput(size_t size, size_t blocks)
U.resize(size * blocks);
V.resize(size * blocks);
_noOfBitSets.resize(blocks);
_noOfOverflowed.resize(blocks);
}
......@@ -421,7 +454,7 @@ GatedFullStokesOutput::GatedFullStokesOutput(size_t nchans, size_t blocks): Outp
G1(nchans, blocks)
{
BOOST_LOG_TRIVIAL(debug) << "Output with " << _blocks << " blocks a " << _nchans << " channels";
_host_power.resize( 8 * ( _nchans * sizeof(IntegratedPowerType) + sizeof(size_t) ) * _blocks);
_host_power.resize( 8 * ( _nchans * sizeof(IntegratedPowerType) + 2 * sizeof(size_t) ) * _blocks);
BOOST_LOG_TRIVIAL(debug) << "Allocated " << _host_power.size() << " bytes.";
};
......@@ -439,7 +472,7 @@ void GatedFullStokesOutput::data2Host(cudaStream_t &_d2h_stream)
for (size_t i = 0; i < G0._noOfBitSets.size(); i++)
{
size_t memslicesize = (_nchans * sizeof(IntegratedPowerType));
size_t memOffset = 8 * i * (memslicesize + sizeof(size_t));
size_t memOffset = 8 * i * (memslicesize + 2 * sizeof(size_t));
// Copy II QQ UU VV
CUDA_ERROR_CHECK(
cudaMemcpyAsync(static_cast<void *>(_host_power.a_ptr() + memOffset) ,
......@@ -495,41 +528,64 @@ for (size_t i = 0; i < G0._noOfBitSets.size(); i++)
static_cast<void *>(G0._noOfBitSets.b_ptr() + i ),
1 * sizeof(size_t),
cudaMemcpyDeviceToHost, _d2h_stream));
std::memcpy(static_cast<void *>(_host_power.a_ptr() + memOffset + 8 * memslicesize + 1 * sizeof(size_t) ),
static_cast<void *>(G0._noOfOverflowed.b_ptr() + i ), 1 * sizeof(size_t) );
CUDA_ERROR_CHECK(
cudaMemcpyAsync( static_cast<void *>(_host_power.a_ptr() + memOffset + 8 * memslicesize + 1 * sizeof(size_t)),
cudaMemcpyAsync( static_cast<void *>(_host_power.a_ptr() + memOffset + 8 * memslicesize + 2 * sizeof(size_t)),
static_cast<void *>(G1._noOfBitSets.b_ptr() + i ),
1 * sizeof(size_t),
cudaMemcpyDeviceToHost, _d2h_stream));
std::memcpy(static_cast<void *>(_host_power.a_ptr() + memOffset + 8 * memslicesize + 3 * sizeof(size_t) ),
static_cast<void *>(G1._noOfOverflowed.b_ptr() + i ), 1 * sizeof(size_t) );
CUDA_ERROR_CHECK(
cudaMemcpyAsync( static_cast<void *>(_host_power.a_ptr() + memOffset + 8 * memslicesize + 2 * sizeof(size_t)),
cudaMemcpyAsync( static_cast<void *>(_host_power.a_ptr() + memOffset + 8 * memslicesize + 4 * sizeof(size_t)),
static_cast<void *>(G0._noOfBitSets.b_ptr() + i ),
1 * sizeof(size_t),
cudaMemcpyDeviceToHost, _d2h_stream));
std::memcpy(static_cast<void *>(_host_power.a_ptr() + memOffset + 8 * memslicesize + 5 * sizeof(size_t) ),
static_cast<void *>(G0._noOfOverflowed.b_ptr() + i ), 1 * sizeof(size_t) );
CUDA_ERROR_CHECK(
cudaMemcpyAsync( static_cast<void *>(_host_power.a_ptr() + memOffset + 8 * memslicesize + 3 * sizeof(size_t)),
cudaMemcpyAsync( static_cast<void *>(_host_power.a_ptr() + memOffset + 8 * memslicesize + 6 * sizeof(size_t)),
static_cast<void *>(G1._noOfBitSets.b_ptr() + i ),
1 * sizeof(size_t),
cudaMemcpyDeviceToHost, _d2h_stream));
std::memcpy(static_cast<void *>(_host_power.a_ptr() + memOffset + 8 * memslicesize + 7 * sizeof(size_t) ),
static_cast<void *>(G1._noOfOverflowed.b_ptr() + i ), 1 * sizeof(size_t) );
CUDA_ERROR_CHECK(
cudaMemcpyAsync( static_cast<void *>(_host_power.a_ptr() + memOffset + 8 * memslicesize + 4 * sizeof(size_t)),
cudaMemcpyAsync( static_cast<void *>(_host_power.a_ptr() + memOffset + 8 * memslicesize + 8 * sizeof(size_t)),
static_cast<void *>(G0._noOfBitSets.b_ptr() + i ),
1 * sizeof(size_t),
cudaMemcpyDeviceToHost, _d2h_stream));
std::memcpy(static_cast<void *>(_host_power.a_ptr() + memOffset + 8 * memslicesize + 9 * sizeof(size_t) ),
static_cast<void *>(G0._noOfOverflowed.b_ptr() + i ), 1 * sizeof(size_t) );
CUDA_ERROR_CHECK(
cudaMemcpyAsync( static_cast<void *>(_host_power.a_ptr() + memOffset + 8 * memslicesize + 5 * sizeof(size_t)),
cudaMemcpyAsync( static_cast<void *>(_host_power.a_ptr() + memOffset + 8 * memslicesize + 10 * sizeof(size_t)),
static_cast<void *>(G1._noOfBitSets.b_ptr() + i ),
1 * sizeof(size_t),
cudaMemcpyDeviceToHost, _d2h_stream));
std::memcpy(static_cast<void *>(_host_power.a_ptr() + memOffset + 8 * memslicesize + 11 * sizeof(size_t) ),
static_cast<void *>(G1._noOfOverflowed.b_ptr() + i ), 1 * sizeof(size_t) );
CUDA_ERROR_CHECK(
cudaMemcpyAsync( static_cast<void *>(_host_power.a_ptr() + memOffset + 8 * memslicesize + 6 * sizeof(size_t)),
cudaMemcpyAsync( static_cast<void *>(_host_power.a_ptr() + memOffset + 8 * memslicesize + 12 * sizeof(size_t)),
static_cast<void *>(G0._noOfBitSets.b_ptr() + i ),
1 * sizeof(size_t),
cudaMemcpyDeviceToHost, _d2h_stream));
std::memcpy(static_cast<void *>(_host_power.a_ptr() + memOffset + 8 * memslicesize + 13 * sizeof(size_t) ),
static_cast<void *>(G0._noOfOverflowed.b_ptr() + i ), 1 * sizeof(size_t) );
CUDA_ERROR_CHECK(
cudaMemcpyAsync( static_cast<void *>(_host_power.a_ptr() + memOffset + 8 * memslicesize + 7 * sizeof(size_t)),
cudaMemcpyAsync( static_cast<void *>(_host_power.a_ptr() + memOffset + 8 * memslicesize + 14 * sizeof(size_t)),
static_cast<void *>(G1._noOfBitSets.b_ptr() + i ),
1 * sizeof(size_t),
cudaMemcpyDeviceToHost, _d2h_stream));
std::memcpy(static_cast<void *>(_host_power.a_ptr() + memOffset + 8 * memslicesize + 15 * sizeof(size_t) ),
static_cast<void *>(G1._noOfOverflowed.b_ptr() + i ), 1 * sizeof(size_t) );
}
}
......
......@@ -276,34 +276,34 @@ TEST(GatedSpectrometer, array_sum) {
}
class GatedTestSink{
class GatedTestSinkSinglePol{
private:
size_t fft_length;
size_t call_count;
size_t nHeaps;
size_t naccumulate;
public:
GatedTestSink(size_t fft_length, size_t nHeaps, size_t naccumulate): fft_length(fft_length), nHeaps(nHeaps), naccumulate(naccumulate), call_count(0) {};
GatedTestSinkSinglePol(size_t fft_length, size_t nHeaps, size_t naccumulate): fft_length(fft_length), nHeaps(nHeaps), naccumulate(naccumulate), call_count(0) {};
// Test the correctness of output of the processing test
void init(psrdada_cpp::RawBytes&){
};
bool operator()(psrdada_cpp::RawBytes& buf)
{
const size_t number_of_sc_items = 1;
const size_t number_of_sc_items = 2;
const size_t nchans = fft_length / 2 + 1;
EXPECT_EQ(buf.used_bytes(), (32 / 8 * nchans + number_of_sc_items* 64 / 8) * 2);
// the first elements of the spectraa are related to the power
float *G0 = reinterpret_cast<float*>(buf.ptr());
float *G1 = reinterpret_cast<float*>(buf.ptr() + nchans * 32 /8);
uint64_t *S = reinterpret_cast<uint64_t*>(buf.ptr() + 2 *nchans * 32 /8);
// Expected half number of samples per gate
uint64_t *S = reinterpret_cast<uint64_t*>(buf.ptr() + 2 *nchans * 32 /8);
EXPECT_EQ(S[0], nHeaps * 4096 / 2);
EXPECT_EQ(S[1], nHeaps * 4096 / 2);
EXPECT_EQ(S[2], nHeaps * 4096 / 2);
//EXPECT_FLOAT_EQ(G0[0], 12.) << "Call count: " << call_count ;
//EXPECT_FLOAT_EQ(G1[1], 23.) << "Call count: " << call_count ;;
// Correct number of overflowed samples
EXPECT_EQ(S[1], 27);
EXPECT_EQ(S[3], 23); // First heap has 23 and bit set, thus G1
call_count ++;
return false;
......@@ -311,7 +311,7 @@ class GatedTestSink{
};
TEST(GatedSpectrometer, processing)
TEST(GatedSpectrometer, processingSinglePol)
{
const size_t nbits = 8;
const size_t nHeaps = 1024;
......@@ -326,9 +326,9 @@ TEST(GatedSpectrometer, processing)
psrdada_cpp::effelsberg::edd::DadaBufferLayout bufferLayout(idbuffer.key(), heapSize, 1);
GatedTestSink sink(fft_length, nHeaps, naccumulate);
GatedTestSinkSinglePol sink(fft_length, nHeaps, naccumulate);
psrdada_cpp::effelsberg::edd::GatedSpectrometer<
GatedTestSink,
GatedTestSinkSinglePol,
psrdada_cpp::effelsberg::edd::SinglePolarizationInput,
psrdada_cpp::effelsberg::edd::GatedPowerSpectrumOutput>
spectrometer(bufferLayout,
......@@ -342,35 +342,114 @@ TEST(GatedSpectrometer, processing)
psrdada_cpp::RawBytes buff(raw_buffer, inputBufferSize, inputBufferSize);
EXPECT_NO_THROW(spectrometer.init(buff));
for (int i = 0; i < nHeaps / 2; i++)
//// fill sci data
uint64_t* sc_items = reinterpret_cast<uint64_t*>(raw_buffer + nHeaps * 4096);
for (int i = 0; i < nHeaps; i+=2)
{
raw_buffer[i] = 0;
sc_items[i] = 0;
sc_items[i+1] = 0;
SET_BIT(sc_items[i], 0);
}
sc_items[0] |= (23UL) << 32;
sc_items[1] |= (27UL) << 32;
// fill data alternating 2,3 per heap
for (int i = 0; i < nHeaps; i+=2)
for (int i = 0; i < 12; i++)
{
for (int j = 0; j < 4096; j++)
{
raw_buffer[i * 4096 + j] = 2;
}
for (int j = 4096; j < 2 * 4096; j++)
{
raw_buffer[i * 4096 + j] = 3;
}
EXPECT_NO_THROW(spectrometer(buff));
}
delete [] raw_buffer;
}
class GatedTestSinkFullStokes{
private:
size_t fft_length;
size_t call_count;
size_t nHeaps;
size_t naccumulate;
public:
GatedTestSinkFullStokes(size_t fft_length, size_t nHeaps, size_t naccumulate): fft_length(fft_length), nHeaps(nHeaps), naccumulate(naccumulate), call_count(0) {};
// Test the correctness of output of the processing test
void init(psrdada_cpp::RawBytes&){
};
bool operator()(psrdada_cpp::RawBytes& buf)
{
const size_t number_of_sc_items = 2;
const size_t nchans = fft_length / 2 + 1;
EXPECT_EQ(buf.used_bytes(), (32 / 8 * nchans + number_of_sc_items* 64 / 8) * 2 * 4); // 4 Spectra (IQUV), 2 Gates
// Expected half number of samples per gate
for(int i =0; i <4; i++)
{
uint64_t *S = reinterpret_cast<uint64_t*>(buf.ptr() + 8 *nchans * 32 /8 + i * 2 *sizeof(size_t));
EXPECT_EQ(S[0], nHeaps * 4096); // expect samples from two polarizations
EXPECT_EQ(S[2], nHeaps * 4096);
// Correct number of overflowed samples
EXPECT_EQ(S[1], 27 + 7);
EXPECT_EQ(S[3], 23 + 3); // First heap has 23 and bit set, thus G1
}
call_count ++;
return false;
};
};
TEST(GatedSpectrometer, processingFullStokes)
{
const size_t nbits = 8;
const size_t nHeaps = 1024;
const size_t fft_length = 1024 * 64;
const size_t naccumulate = 4096 * nHeaps / fft_length;
const size_t heapSize = 4096 * nbits / 8;
const size_t inputBufferSize = 2 * nHeaps * (heapSize + 64 / 8) ;
psrdada_cpp::DadaDB idbuffer(5, inputBufferSize, 1, 4096);
idbuffer.create();
psrdada_cpp::effelsberg::edd::DadaBufferLayout bufferLayout(idbuffer.key(), heapSize, 1);
GatedTestSinkFullStokes sink(fft_length, nHeaps, naccumulate);
psrdada_cpp::effelsberg::edd::GatedSpectrometer<
GatedTestSinkFullStokes,
psrdada_cpp::effelsberg::edd::DualPolarizationInput,
psrdada_cpp::effelsberg::edd::GatedFullStokesOutput>
spectrometer(bufferLayout,
0, 0,
fft_length,
naccumulate, nbits,
sink);
char *raw_buffer = new char[inputBufferSize];
psrdada_cpp::RawBytes buff(raw_buffer, inputBufferSize, inputBufferSize);
EXPECT_NO_THROW(spectrometer.init(buff));
//// fill sci data
uint64_t* sc_items = reinterpret_cast<uint64_t*>(raw_buffer + nHeaps * 4096);
for (int i = 0; i < nHeaps; i+=2)
uint64_t* sc_items = reinterpret_cast<uint64_t*>(raw_buffer + 2*nHeaps * 4096);
for (int i = 0; i < 2 * nHeaps; i+=4)
{
sc_items[i] = 0;
sc_items[i+1] = 0;
sc_items[i+2] = 0;
sc_items[i+3] = 0;
SET_BIT(sc_items[i], 0);
//SET_BIT(sc_items[i+1], 0);
SET_BIT(sc_items[i + 1], 0);
}
sc_items[0] |= (23UL) << 32;
sc_items[1] |= (3UL) << 32;
sc_items[2] |= (27UL) << 32;
sc_items[3] |= (7UL) << 32;