Merge branch 'develop' of ssh://code.ascee.nl:12001/ASCEE/lasp into develop

This commit is contained in:
Casper Jansen 2022-10-13 12:07:05 +02:00
commit a73863d044
8 changed files with 222 additions and 131 deletions

View File

@ -8,9 +8,10 @@ set(CMAKE_CXX_STANDARD_REQUIRED)
option(LASP_DOUBLE_PRECISION "Compile as double precision floating point" ON) option(LASP_DOUBLE_PRECISION "Compile as double precision floating point" ON)
option(LASP_HAS_RTAUDIO "Compile with RtAudio Daq backend" ON) option(LASP_HAS_RTAUDIO "Compile with RtAudio Daq backend" ON)
option(LASP_HAS_ULDAQ "Compile with UlDaq backend" ON) option(LASP_HAS_ULDAQ "Compile with UlDaq backend" ON)
option(LASP_BUILD_TUNED "Tune build for current machine" OFF) option(LASP_BUILD_TUNED "Tune build for current machine (Experimental / untested)" OFF)
option(LASP_WITH_OPENMP "Use OpenMP parallelization" ON) option(LASP_WITH_OPENMP "Use OpenMP parallelization (Experimental: crashes SHOULD BE EXPECTED)" OFF)
set(LASP_MAX_NFFT "33554432" CACHE STRING "Max FFT size") set(LASP_MAX_NFFT "33554432" CACHE STRING "Max FFT size")
option(LASP_BUILD_CPP_TESTS "Build CPP test code" OFF)
# Use ccache if available # Use ccache if available
find_program(CCACHE_PROGRAM ccache) find_program(CCACHE_PROGRAM ccache)
@ -100,4 +101,7 @@ include(rtaudio)
include(uldaq) include(uldaq)
# #
add_subdirectory(src/lasp) add_subdirectory(src/lasp)
if(LASP_BUILD_CPP_TESTS)
add_subdirectory(test)
endif()

View File

@ -6,7 +6,7 @@ add_definitions(-DARMA_DONT_USE_WRAPPER)
configure_file(lasp_config.h.in lasp_config.h) configure_file(lasp_config.h.in lasp_config.h)
include_directories(${CMAKE_CURRENT_BINARY_DIR}) include_directories(${CMAKE_CURRENT_BINARY_DIR})
include_directories(SYSTEM include_directories(SYSTEM
../../third_party/armadillo-code/include) ${PROJECT_SOURCE_DIR}/third_party/armadillo-code/include)
include_directories(../../third_party/DebugTrace-cpp/include) include_directories(../../third_party/DebugTrace-cpp/include)
include_directories(../../third_party/gsl-lite/include) include_directories(../../third_party/gsl-lite/include)

View File

@ -64,17 +64,18 @@ void DaqData::copyInFromRaw(const std::vector<byte_t *> &ptrs) {
/* std::copy(ptr, ptr + sw * nframes, &_data[sw * ch * nframes]); */ /* std::copy(ptr, ptr + sw * nframes, &_data[sw * ch * nframes]); */
ch++; ch++;
} }
dmat data2(nframes, nchannels);
} }
void DaqData::copyToRaw(const us channel, byte_t *ptr) { void DaqData::copyToRaw(const us channel, byte_t *ptr) {
/* std::copy(raw_ptr(0, channel), raw_ptr(nframes, channel), ptr); */ /* std::copy(raw_ptr(0, channel), raw_ptr(nframes, channel), ptr); */
assert(channel<nchannels);
assert(ptr);
memcpy(ptr, raw_ptr(0, channel), sw*nframes); memcpy(ptr, raw_ptr(0, channel), sw*nframes);
} }
template <typename T> template <typename T>
d DaqData::toFloat(const us frame, const us channel) const { d DaqData::toFloat(const us frame, const us channel) const {
DEBUGTRACE_ENTER; /* DEBUGTRACE_ENTER; */
if constexpr (std::is_integral<T>::value) { if constexpr (std::is_integral<T>::value) {
return static_cast<d>(value<T>(frame, channel)) / return static_cast<d>(value<T>(frame, channel)) /
std::numeric_limits<T>::max(); std::numeric_limits<T>::max();

View File

@ -1,34 +1,30 @@
// lasp_fft.c /* #define DEBUGTRACE_ENABLED */
//
// Author: J.A. de Jong - ASCEE
//
// Description:
// FFt implementation
//
//////////////////////////////////////////////////////////////////////
#include <memory> #include <memory>
#include <cassert> #include <cassert>
/* #define DEBUGTRACE_ENABLED */
#include "lasp_fft.h"
#include "debugtrace.hpp" #include "debugtrace.hpp"
#include "lasp_config.h" #include "lasp_config.h"
#include "lasp_fft.h"
using rte = std::runtime_error; using rte = std::runtime_error;
#if LASP_FFT_BACKEND == Armadillo #if LASP_FFT_BACKEND == Armadillo
#include "fftpack.h"
class Fft_impl { class Fft_impl {
public: public:
us nfft; us nfft;
Fft_impl(const us nfft) : nfft(nfft) { Fft_impl(const us nfft) : nfft(nfft) { DEBUGTRACE_ENTER; }
throw runtime_error(
"This code does not output correct results, as it computes the full "
"FFT. It needs to be reworked to single-sided spectra.");
DEBUGTRACE_ENTER; vc fft(const vd &time) {
} DEBUGTRACE_ENTER;
vc res = arma::fft(time);
vd ifft(const vc &f) { return arma::ifft(f); } return res.rows(0, nfft / 2);
vc fft(const vd &time) { return arma::fft(time); } }
vd ifft(const vc &f) {
DEBUGTRACE_ENTER;
vc f_full(nfft);
vc f_above_nyq = arma::conj(arma::reverse(f));
f_full.rows(0, nfft / 2) = f;
f_full.rows(nfft / 2, nfft - 1) = f_above_nyq.rows(0, nfft/2-1);
return arma::real(arma::ifft(f_full));
}
}; };
#elif LASP_FFT_BACKEND == FFTW #elif LASP_FFT_BACKEND == FFTW
#include <fftw3.h> #include <fftw3.h>
@ -37,77 +33,76 @@ class Fft_impl {
* @brief FFTW implementation * @brief FFTW implementation
*/ */
class Fft_impl { class Fft_impl {
public: public:
us nfft; us nfft;
fftw_plan forward_plan = nullptr; fftw_plan forward_plan = nullptr;
fftw_plan reverse_plan = nullptr; fftw_plan reverse_plan = nullptr;
fftw_complex *frequencyDomain = nullptr; fftw_complex *frequencyDomain = nullptr;
d *timeDomain = nullptr; d *timeDomain = nullptr;
Fft_impl(const us nfft) : nfft(nfft) { Fft_impl(const us nfft) : nfft(nfft) {
timeDomain = (d *)fftw_malloc(sizeof(d) * nfft); timeDomain = (d *)fftw_malloc(sizeof(d) * nfft);
frequencyDomain = frequencyDomain =
(fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (nfft / 2 + 1)); (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (nfft / 2 + 1));
forward_plan = forward_plan =
fftw_plan_dft_r2c_1d(nfft, timeDomain, frequencyDomain, FFTW_MEASURE); fftw_plan_dft_r2c_1d(nfft, timeDomain, frequencyDomain, FFTW_MEASURE);
reverse_plan = reverse_plan =
fftw_plan_dft_c2r_1d(nfft, frequencyDomain, timeDomain, FFTW_MEASURE); fftw_plan_dft_c2r_1d(nfft, frequencyDomain, timeDomain, FFTW_MEASURE);
if(!forward_plan || !reverse_plan || !timeDomain || !frequencyDomain) { if (!forward_plan || !reverse_plan || !timeDomain || !frequencyDomain) {
throw rte("Error allocating FFT"); throw rte("Error allocating FFT");
} }
};
}; ~Fft_impl() {
fftw_destroy_plan(forward_plan);
fftw_destroy_plan(reverse_plan);
fftw_free(frequencyDomain);
fftw_free(timeDomain);
}
~Fft_impl() { vc fft(const vd &time) {
fftw_destroy_plan(forward_plan); if (time.n_rows != nfft) {
fftw_destroy_plan(reverse_plan); throw rte("Invalid size of input vector, should be equal to nfft");
fftw_free(frequencyDomain);
fftw_free(timeDomain);
} }
vc fft(const vd &time) { vc res(nfft / 2 + 1);
if(time.n_rows != nfft) { memcpy(timeDomain, time.memptr(), sizeof(d) * nfft);
throw rte("Invalid size of input vector, should be equal to nfft");
}
vc res(nfft / 2 + 1); fftw_execute(forward_plan);
memcpy(timeDomain, time.memptr(), sizeof(d) * nfft);
fftw_execute(forward_plan); memcpy(reinterpret_cast<void *>(res.memptr()), frequencyDomain,
sizeof(c) * (nfft / 2 + 1));
memcpy(reinterpret_cast<void *>(res.memptr()), frequencyDomain, return res;
sizeof(c) * (nfft / 2 + 1)); }
vd ifft(const vc &f) {
return res; DEBUGTRACE_ENTER;
if (f.n_rows != nfft / 2 + 1) {
throw rte("Invalid size of input vector, should be equal to nfft/2+1");
} }
vd ifft(const vc &f) { memcpy(frequencyDomain,
DEBUGTRACE_ENTER; reinterpret_cast<void *>(const_cast<c *>(f.memptr())),
if(f.n_rows != nfft/2+1) { (nfft / 2 + 1) * sizeof(c));
throw rte("Invalid size of input vector, should be equal to nfft/2+1");
}
memcpy(frequencyDomain,
reinterpret_cast<void *>(const_cast<c *>(f.memptr())),
(nfft / 2 + 1) * sizeof(c));
fftw_execute(reverse_plan); fftw_execute(reverse_plan);
vd res(nfft); vd res(nfft);
memcpy(res.memptr(), timeDomain, nfft * sizeof(d)); memcpy(res.memptr(), timeDomain, nfft * sizeof(d));
/* Scale by dividing by nfft. Checked with numpy implementation /* Scale by dividing by nfft. Checked with numpy implementation
* that this indeed needs to be done for FFTW. */ * that this indeed needs to be done for FFTW. */
res *= 1.0 / nfft; res *= 1.0 / nfft;
return res; return res;
} }
}; };
#else #else
#error \ #error \
"Cannot compile lasp_ffc.c, no FFT backend specified. Should either be FFTPack, or FFTW" "Cannot compile lasp_ffc.c, no FFT backend specified. Should either be FFTPack, or FFTW"
#endif #endif
Fft::Fft(const us nfft) { Fft::Fft(const us nfft) {
@ -117,20 +112,26 @@ Fft::Fft(const us nfft) {
} }
if (nfft >= LASP_MAX_NFFT) { if (nfft >= LASP_MAX_NFFT) {
throw rte("Invalid nfft, should be smaller than: " + throw rte("Invalid nfft, should be smaller than: " +
std::to_string(LASP_MAX_NFFT)); std::to_string(LASP_MAX_NFFT));
} }
_impl = std::make_unique<Fft_impl>(nfft); _impl = std::make_unique<Fft_impl>(nfft);
} }
Fft::~Fft() { } Fft::~Fft() {}
us Fft::nfft() const { assert(_impl); return _impl->nfft; } us Fft::nfft() const {
assert(_impl);
return _impl->nfft;
}
vc Fft::fft(const vd &timedata) { assert(_impl); return _impl->fft(timedata); } vc Fft::fft(const vd &timedata) {
assert(_impl);
return _impl->fft(timedata);
}
cmat Fft::fft(const dmat &freqdata) { cmat Fft::fft(const dmat &freqdata) {
DEBUGTRACE_ENTER; DEBUGTRACE_ENTER;
assert(_impl); assert(_impl);
cmat res(_impl->nfft/2+1, freqdata.n_cols); cmat res(_impl->nfft / 2 + 1, freqdata.n_cols);
/// * WARNING *. This was source of a serious bug. It is not possible to run /// * WARNING *. This was source of a serious bug. It is not possible to run
/// FFT's and IFFT's on the same _impl, as it overwrites the same memory. /// FFT's and IFFT's on the same _impl, as it overwrites the same memory.
/// Uncommenting the line below results in faulty results. /// Uncommenting the line below results in faulty results.
@ -143,7 +144,7 @@ cmat Fft::fft(const dmat &freqdata) {
vd Fft::ifft(const vc &freqdata) { vd Fft::ifft(const vc &freqdata) {
DEBUGTRACE_ENTER; DEBUGTRACE_ENTER;
assert(_impl); assert(_impl);
return _impl->ifft(freqdata); return _impl->ifft(freqdata);
} }
dmat Fft::ifft(const cmat &freqdata) { dmat Fft::ifft(const cmat &freqdata) {
@ -159,7 +160,7 @@ dmat Fft::ifft(const cmat &freqdata) {
return res; return res;
} }
void Fft::load_fft_wisdom(const std::string& wisdom) { void Fft::load_fft_wisdom(const std::string &wisdom) {
#if LASP_FFT_BACKEND == Armadillo #if LASP_FFT_BACKEND == Armadillo
#elif LASP_FFT_BACKEND == FFTW #elif LASP_FFT_BACKEND == FFTW
if (wisdom.length() > 0) { if (wisdom.length() > 0) {
@ -179,8 +180,8 @@ std::string Fft::store_fft_wisdom() {
// put it in this container by copying in. Not a good solution if this has to // put it in this container by copying in. Not a good solution if this has to
// happen often. // happen often.
// Fortunately, this function is only called at the end of the program. // Fortunately, this function is only called at the end of the program.
char* wis = fftw_export_wisdom_to_string(); char *wis = fftw_export_wisdom_to_string();
std::string res {wis}; std::string res{wis};
free(wis); free(wis);
return res; return res;
#endif #endif

View File

@ -1,4 +1,4 @@
#define DEBUGTRACE_ENABLED /* #define DEBUGTRACE_ENABLED */
#include "lasp_rtaps.h" #include "lasp_rtaps.h"
#include "debugtrace.hpp" #include "debugtrace.hpp"
#include <mutex> #include <mutex>
@ -13,9 +13,9 @@ RtAps::RtAps(StreamMgr &mgr, const Filter *freqWeightingFilter,
: ThreadedInDataHandler(mgr), : ThreadedInDataHandler(mgr),
_ps(nfft, w, overlap_percentage, time_constant) { _ps(nfft, w, overlap_percentage, time_constant) {
/* if (freqWeightingFilter != nullptr) { */ if (freqWeightingFilter != nullptr) {
/* _filterPrototype = freqWeightingFilter->clone(); */ _filterPrototype = freqWeightingFilter->clone();
/* } */ }
start(); start();
} }
@ -29,30 +29,29 @@ bool RtAps::inCallback_threaded(const DaqData &data) {
std::scoped_lock<std::mutex> lck(_mtx); std::scoped_lock<std::mutex> lck(_mtx);
dmat fltdata = data.toFloat(); dmat fltdata = data.toFloat();
data.print(); const us nchannels = fltdata.n_cols;
/* const us nchannels = fltdata.n_cols; */
/* if (_filterPrototype) { */ if (_filterPrototype) {
/* // Adjust number of filters, if necessary */ // Adjust number of filters, if necessary
/* if (nchannels > _freqWeightingFilter.size()) { */ if (nchannels > _freqWeightingFilter.size()) {
/* while (nchannels > _freqWeightingFilter.size()) { */ while (nchannels > _freqWeightingFilter.size()) {
/* _freqWeightingFilter.emplace_back(_filterPrototype->clone()); */ _freqWeightingFilter.emplace_back(_filterPrototype->clone());
/* } */ }
/* for (auto &filter : _freqWeightingFilter) { */ for (auto &filter : _freqWeightingFilter) {
/* filter->reset(); */ filter->reset();
/* } */ }
/* } */ }
/* // Apply filtering */ // Apply filtering
/* #pragma omp parallel for */ #pragma omp parallel for
/* for (us i = 0; i < nchannels; i++) { */ for (us i = 0; i < nchannels; i++) {
/* vd col = fltdata.col(i); */ vd col = fltdata.col(i);
/* _freqWeightingFilter.at(i)->filter(col); */ _freqWeightingFilter.at(i)->filter(col);
/* fltdata.col(i) = col; */ fltdata.col(i) = col;
/* } */ }
/* } // End of if(_filterPrototype) */ } // End of if(_filterPrototype)
_ps.compute(fltdata); _ps.compute(fltdata);
@ -75,9 +74,4 @@ std::unique_ptr<ccube> RtAps::getCurrentValue() {
auto est = _ps.get_est(); auto est = _ps.get_est();
return std::make_unique<ccube>(est.value_or(ccube())); return std::make_unique<ccube>(est.value_or(ccube()));
/* return std::move(_latest_est); */
/* if (_latest_est) { */
/* return std::make_unique<cube>(cube(*_latest_est)); */
/* } else */
/* return nullptr; */
} }

View File

@ -1,17 +1,8 @@
include_directories(${CMAKE_SOURCE_DIR}/lasp/c) include_directories(SYSTEM
include_directories(${CMAKE_SOURCE_DIR}/lasp/device) ${PROJECT_SOURCE_DIR}/third_party/armadillo-code/include)
add_executable(test_daq test_daq.cpp)
# add_executable(test_bf test_bf.c) add_executable(test_smgr test_smgr.cpp)
# add_executable(test_workers test_workers.c) include_directories(../src/lasp/device ../src/lasp/dsp ../src/lasp)
# add_executable(test_fft test_fft.c) include_directories(../third_party/gsl-lite/include)
# add_executable(test_math test_math.c) target_link_libraries(test_daq lasp_device_lib lasp_dsp_lib)
target_link_libraries(test_smgr lasp_device_lib lasp_dsp_lib)
# target_link_libraries(test_bf lasp_lib ${LASP_THREADING_LIBRARIES})
# target_link_libraries(test_fft lasp_lib ${LASP_THREADING_LIBRARIES})
# target_link_libraries(test_workers lasp_lib ${LASP_THREADING_LIBRARIES})
# target_link_libraries(test_math lasp_lib ${LASP_THREADING_LIBRARIES})
if(LASP_ULDAQ)
add_executable(test_uldaq test_uldaq.cpp)
target_link_libraries(test_uldaq cpp_daq)
endif(LASP_ULDAQ)

48
test/test_daq.cpp Normal file
View File

@ -0,0 +1,48 @@
#include "lasp_daq.h"
#include "lasp_deviceinfo.h"
#include <memory>
#include <vector>
using namespace std;
bool inCallback(const DaqData& d) {
d.toFloat();
return true;
}
int main(int argc, const char **const argv) {
std::vector<DeviceInfo> devs = DeviceInfo::getDeviceInfo();
DeviceInfo *mon_device = nullptr;
for (auto &d : devs) {
string name_lower = d.device_name;
transform(name_lower.begin(), name_lower.end(), name_lower.begin(),
::tolower);
if (name_lower.find("monitor") != string::npos) {
mon_device = &d;
}
}
if (!mon_device) {
cerr << "Could not find monitor device\n";
exit(1);
}
DaqConfiguration config(*mon_device);
config.inchannel_config[0].enabled = true;
config.inchannel_config[1].enabled = true;
unique_ptr<Daq> daq = Daq::createDaq(*mon_device, config);
InDaqCallback icb = inCallback;
OutDaqCallback ocb;
daq->start(icb, ocb);
cout << "Press <enter> to stop" << endl;
cin.get();
return 0;
}

52
test/test_smgr.cpp Normal file
View File

@ -0,0 +1,52 @@
#include "lasp_daq.h"
#include "lasp_deviceinfo.h"
#include "lasp_streammgr.h"
#include "lasp_ppm.h"
#include <memory>
#include <vector>
using namespace std;
bool inCallback(const DaqData& d) {
d.toFloat();
return true;
}
int main(int argc, const char **const argv) {
StreamMgr& mgr = StreamMgr::getInstance();
std::vector<DeviceInfo> devs = mgr.getDeviceInfo();
DeviceInfo *mon_device = nullptr;
for (auto &d : devs) {
string name_lower = d.device_name;
transform(name_lower.begin(), name_lower.end(), name_lower.begin(),
::tolower);
if (name_lower.find("monitor") != string::npos) {
mon_device = &d;
}
}
if (!mon_device) {
cerr << "Could not find monitor device\n";
exit(1);
}
cout << "Found device. Name: " << mon_device->device_name << endl;
DaqConfiguration config(*mon_device);
config.inchannel_config.at(0).enabled = true;
config.inchannel_config.at(1).enabled = true;
mgr.startStream(config);
PPMHandler ppm(mgr);
ppm.start();
cout << "Press <enter> to stop" << endl;
cin.get();
return 0;
}