ASCEE/lasp
ASCEE/lasp
5
1
Fork 0
Code Issues 4 Pull Requests Actions Releases 12 Activity

Sweep bug fixed. There is still something weird with RtAudio: only one channel is outputting signal.

Browse Source
This commit is contained in:
Anne de Jong 2022-09-22 21:02:41 +02:00
parent b200b465f6
commit 5ce5fba50b
12 changed files with 228 additions and 83 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
CMakeLists.txt
View File

@ -9,6 +9,7 @@ option(LASP_DOUBLE_PRECISION "Compile as double precision floating point" ON)
option(LASP_HAS_RTAUDIO "Compile with RtAudio Daq backend" ON)
option(LASP_HAS_ULDAQ "Compile with UlDaq backend" ON)
option(LASP_BUILD_TUNED "Tune build for current machine" OFF)
option(LASP_WITH_OPENMP "Use OpenMP parallelization" ON)
# Use ccache if available
find_program(CCACHE_PROGRAM ccache)
@ -42,6 +43,8 @@ set(CMAKE_POSITION_INDEPENDENT_CODE ON)
if(${CMAKE_BUILD_TYPE} STREQUAL "Debug")
set(LASP_DEBUG True)
# This does not work nicely with RtAudio
# add_definitions(-D_GLIBCXX_DEBUG)
else()
set(LASP_DEBUG False)
endif()
@ -65,8 +68,12 @@ add_definitions(-DLASP_MAX_NFFT=33554432) # 2**25
# ####################################### End of user-adjustable variables section
include(OSSpecific)
# Enable openmp
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS} -Wall -Wextra -Wno-type-limits")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -Wno-type-limits")
# Enable openmp, if set
if(LASP_WITH_OPENMP)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
endif()
set(CMAKE_C_FLAGS_RELEASE "-O3 -flto -mfpmath=sse -march=x86-64 -mtune=native \
-fdata-sections -ffunction-sections -fomit-frame-pointer -finline-functions")
@ -78,7 +85,6 @@ add_subdirectory(third_party/carma)
if(LASP_FFT_BACKEND STREQUAL "FFTW")
find_library(fftw3 REQUIRED NAMES fftw fftw3)
set(LASP_FFT_LIBS "fftw3")
message("LASP_FFT_LIBS = ${LASP_FFT_LIBS}")
elseif(LASP_FFT_BACKEND STREQUAL "Armadillo")
endif()

1
src/lasp/CMakeLists.txt
View File

@ -28,6 +28,7 @@ pybind11_add_module(lasp_cpp MODULE lasp_cpp.cpp
pybind11/lasp_daq.cpp
pybind11/lasp_deviceinfo.cpp
pybind11/lasp_pyindatahandler.cpp
pybind11/lasp_siggen.cpp
)
target_link_libraries(lasp_cpp PRIVATE lasp_device_lib lasp_dsp_lib

8
src/lasp/device/lasp_daqdata.h
View File

@ -60,7 +60,7 @@ public:
*/
DaqData() : DaqData(0, 0, DataTypeDescriptor::DataType::dtype_int8) {}
virtual ~DaqData() = default;
DaqData& operator=(const DaqData&) = default;
DaqData &operator=(const DaqData &) = default;
/**
* @brief Return pointer to the raw data corresponding to a certain sample.
@ -74,6 +74,12 @@ public:
return &(_data.data()[sw * (frame + channel * nframes)]);
}
void setSlow(const us frame, const us channel, const int8_t* val) {
for(us i=0;i<sw;i++) {
_data.at(i+sw*(frame + channel*nframes)) = val[i];
}
}
/**
* @brief Return the total number of bytes
*

8
src/lasp/device/lasp_rtaudiodaq.cpp
View File

@ -327,14 +327,14 @@ class RtAudioDaq : public Daq {
if (outputBuffer) {
std::vector<uint8_t *> ptrs;
ptrs.reserve(neninchannels);
DaqData data(nenoutchannels, nFramesPerBlock, dtype);
ptrs.reserve(nenoutchannels);
/* outCallback */
for (int ch = 0; ch <= getHighestOutChannel(); ch++) {
/* DEBUGTRACE_PRINT(outchannel_config.at(ch).enabled); */
if (outchannel_config.at(ch).enabled) {
ptrs.push_back(&static_cast<uint8_t *>(
outputBuffer)[sw * nenoutchannels * ch * nFramesPerBlock]);
ptrs.push_back(&(static_cast<uint8_t *>(
outputBuffer)[sw * nenoutchannels * ch * nFramesPerBlock]));
}
}
DaqData d{nenoutchannels, nFramesPerBlock, dtype};

61
src/lasp/device/lasp_streammgr.cpp
View File

@ -89,6 +89,25 @@ bool StreamMgr::inCallback(const DaqData &data) {
return true;
}
void StreamMgr::setSiggen(std::shared_ptr<Siggen> siggen) {
DEBUGTRACE_ENTER;
checkRightThread();
std::scoped_lock<std::mutex> lck(_siggen_mtx);
// If not set to nullptr, and a stream is running, we update the signal
// generator by resetting it.
if (isStreamRunningOK(StreamType::output) && siggen) {
const Daq *daq = getDaq(StreamType::output);
assert(daq != nullptr);
// Reset the signal generator.
siggen->reset(daq->samplerate());
}
_siggen = siggen;
}
#define DEBUG_FILLDATA 0
/**
* @brief Converts from double precision floating point to output signal in
* non-interleaving format.
@ -100,48 +119,45 @@ bool StreamMgr::inCallback(const DaqData &data) {
* @return
*/
template <typename T> bool fillData(DaqData &data, const vd &signal) {
/* DEBUGTRACE_ENTER; */
assert(data.nframes == signal.size());
/* cerr << "SFSG: data.nframes:" << data.nframes << endl; */
/* cerr << "SFSG: data.nchannels:" << data.nchannels << endl; */
T *res = reinterpret_cast<T *>(data.raw_ptr());
if (std::is_floating_point<T>()) {
for (us ch = 0; ch < data.nchannels; ch++) {
for (us frame = 0; frame < data.nframes; frame++) {
#if DEBUG_FILLDATA == 1
DEBUGTRACE_PRINT("SLOW flt");
data.setSlow(frame, ch,
reinterpret_cast<const int8_t *>(&signal[frame]));
#else
res[ch * data.nframes + frame] = signal[frame];
#endif
}
}
} else {
for (us ch = 0; ch < data.nchannels; ch++) {
for (us frame = 0; frame < data.nframes; frame++) {
res[ch * data.nframes + frame] =
(signal[frame] * std::numeric_limits<T>::max());
const T val = (signal[frame] * std::numeric_limits<T>::max());
#if DEBUG_FILLDATA == 1
data.setSlow(frame, ch, reinterpret_cast<const int8_t *>(&val));
#else
res[ch * data.nframes + frame] = val;
#endif
}
}
}
return true;
}
void StreamMgr::setSiggen(std::shared_ptr<Siggen> siggen) {
checkRightThread();
std::scoped_lock<std::mutex> lck(_siggen_mtx);
_siggen = siggen;
// If not set to nullptr, and a stream is running, we update the signal
// generator by resetting it.
if (isStreamRunningOK(StreamType::output) && siggen) {
const Daq *daq = getDaq(StreamType::output);
assert(daq != nullptr);
// Reset the signal generator.
_siggen->reset(daq->samplerate());
}
}
bool StreamMgr::outCallback(DaqData &data) {
/* DEBUGTRACE_ENTER; */
std::scoped_lock<std::mutex> lck(_siggen_mtx);
if (_siggen) {
vd signal = _siggen->genSignal(data.nframes);
switch (data.dtype) {
@ -237,6 +253,8 @@ void StreamMgr::startStream(const DaqConfiguration &config) {
if (isOutput) {
if (_siggen) {
DEBUGTRACE_PRINT("Resetting _siggen with new samplerate of ");
DEBUGTRACE_PRINT(daq->samplerate());
_siggen->reset(daq->samplerate());
}
outCallback = std::bind(&StreamMgr::outCallback, this, _1);
@ -248,6 +266,9 @@ void StreamMgr::startStream(const DaqConfiguration &config) {
if (isInput) {
_inputStream = std::move(daq);
if (daq->duplexMode()) {
assert(!_outputStream);
}
} else {
assert(isOutput);
_outputStream = std::move(daq);
@ -291,7 +312,7 @@ void StreamMgr::removeInDataHandler(InDataHandler &handler) {
}
Daq::StreamStatus StreamMgr::getStreamStatus(const StreamType type) const {
DEBUGTRACE_ENTER;
/* DEBUGTRACE_ENTER; */
checkRightThread();
// Default constructor, says stream is not running, but also no errors

24
src/lasp/dsp/lasp_siggen.cpp
View File

@ -1,21 +1,25 @@
/* #define DEBUGTRACE_ENABLED */
#include "debugtrace.hpp"
#include "lasp_siggen.h"
#include "lasp_mathtypes.h"
#include <cassert>
#include <type_traits>
using std::cerr;
using std::endl;
inline d level_amp(d level_dB){
return pow(10, level_dB/20);
}
inline d level_amp(d level_dB) { return pow(10, level_dB / 20); }
using mutexlock = std::scoped_lock<std::mutex>;
vd Siggen::genSignal(const us nframes) {
DEBUGTRACE_ENTER;
mutexlock lck(_mtx);
vd signal(nframes, arma::fill::value(_dc_offset));
if (!_muted) {
vd signal_dynamic = _level_linear*genSignalUnscaled(nframes);
if(_filter) {
vd signal_dynamic = _level_linear * genSignalUnscaled(nframes);
if (_filter) {
_filter->filter(signal_dynamic);
}
signal += signal_dynamic;
@ -23,20 +27,24 @@ vd Siggen::genSignal(const us nframes) {
return signal;
}
void Siggen::setFilter(std::shared_ptr<Filter>& filter) {
void Siggen::setFilter(std::shared_ptr<Filter> filter) {
DEBUGTRACE_ENTER;
mutexlock lck(_mtx);
_filter = filter;
}
void Siggen::setDCOffset(const d offset) {
DEBUGTRACE_ENTER;
mutexlock lck(_mtx);
_dc_offset = offset;
}
void Siggen::setLevel(const d level,bool dB) {
void Siggen::setLevel(const d level, bool dB) {
DEBUGTRACE_ENTER;
mutexlock lck(_mtx);
_level_linear = dB ? level_amp(level) : level;
}
void Siggen::reset(const d newFs) {
DEBUGTRACE_ENTER;
mutexlock lck(_mtx);
fs = newFs;
_fs = newFs;
resetImpl();
}

24
src/lasp/dsp/lasp_siggen.h
View File

@ -8,17 +8,27 @@ class StreamMgr;
class DaqData;
/**
* @brief Signal generation base class
* \addtogroup dsp
* @{
*/
/**
* \defgroup siggen Signal generators
* @{
*/
/**
* @brief Signal generation abstract base class. Implementation is required for
* resetImpl(), genSignalUnscaled() and destructor.
*/
class Siggen {
private:
std::shared_ptr<Filter> _filter;
d _dc_offset = 0, _level_linear = 1;
bool _muted = false;
std::mutex _mtx;
protected:
d fs = -1;
std::mutex _mtx;
d _fs = 0;
virtual void resetImpl() = 0;
virtual vd genSignalUnscaled(const us nframes) = 0;
@ -32,7 +42,7 @@ public:
*
* @param f The filter to install.
*/
void setFilter(std::shared_ptr<Filter>& f);
void setFilter(std::shared_ptr<Filter> f);
/**
* @brief Set a linear DC offset value to the signal
@ -42,7 +52,8 @@ public:
void setDCOffset(d offset);
/**
* @brief Mute the signal. Passes through the DC offset.
* @brief Mute the signal. Passes through the DC offset. No lock is hold. If
* it just works one block later, than that is just the case.
*
* @param mute if tre
*/
@ -75,3 +86,6 @@ public:
*/
vd genSignal(const us nframes);
};
/** @} */
/** @} */

65
src/lasp/dsp/lasp_siggen_impl.cpp
View File

@ -5,11 +5,14 @@
// Description:
// Signal generators implementation
//////////////////////////////////////////////////////////////////////
/* #define TRACERPLUS (-5) */
#define DEBUGTRACE_ENABLED
#include "debugtrace.hpp"
#include "lasp_siggen_impl.h"
#include "debugtrace.hpp"
#include "lasp_mathtypes.h"
using rte = std::runtime_error;
DEBUGTRACE_VARIABLES;
/** The fixed number of Newton iterations t.b.d. for tuning the sweep start and
@ -23,13 +26,14 @@ vd Noise::genSignalUnscaled(us nframes) {
}
void Noise::resetImpl() {}
Sine::Sine(const d freq) : omg(2 * arma::datum::pi * freq) {}
Sine::Sine(const d freq) : omg(2 * arma::datum::pi * freq) { DEBUGTRACE_ENTER;}
vd Sine::genSignalUnscaled(const us nframes) {
/* DEBUGTRACE_ENTER; */
const d pi = arma::datum::pi;
vd phase_vec =
arma::linspace(phase, phase + omg * (nframes - 1) / fs, nframes);
phase += omg * nframes / fs;
arma::linspace(phase, phase + omg * (nframes - 1) / _fs, nframes);
phase += omg * nframes / _fs;
while (phase > 2 * arma::datum::pi) {
phase -= 2 * pi;
}
@ -37,10 +41,13 @@ vd Sine::genSignalUnscaled(const us nframes) {
}
vd Periodic::genSignalUnscaled(const us nframes) {
us signal_idx = 0;
vd res(nframes);
while (signal_idx < nframes) {
res[signal_idx] = _signal[_cur_pos];
if(_signal.size() == 0) {
throw rte("No signal defined while calling");
}
for(us i=0;i<nframes;i++) {
res(i) = _signal[_cur_pos];
_cur_pos++;
_cur_pos %= _signal.size();
}
@ -48,29 +55,35 @@ vd Periodic::genSignalUnscaled(const us nframes) {
}
Sweep::Sweep(const d fl, const d fu, const d Ts, const d Tq, const us flags)
: fl_(fl), fu_(fu), Ts(Ts), Tq(Tq), flags(flags) {
if (fl <= 0 || fu < fl || Ts <= 0) {
throw std::runtime_error("Invalid sweep parameters");
: fl_(fl), fu_(fu), Ts(Ts), Tq(Tq), flags(flags) {
if (fl <= 0 || fu < fl || Ts <= 0) {
throw std::runtime_error("Invalid sweep parameters");
}
if ((flags & ForwardSweep) && (flags & BackwardSweep)) {
throw std::runtime_error(
"Both forward and backward sweep flag set. Please only set either one "
"or none for a continuous sweep");
}
if ((flags & LinearSweep) && (flags & LogSweep)) {
throw std::runtime_error(
"Both logsweep and linear sweep flag set. Please only set either one.");
}
}
if ((flags & ForwardSweep) && (flags & BackwardSweep)) {
throw std::runtime_error(
"Both forward and backward sweep flag set. Please only set either one "
"or none for a continuous sweep");
}
}
void Sweep::resetImpl() {
DEBUGTRACE_ENTER;
_cur_pos = 0;
bool forward_sweep = flags & ForwardSweep;
bool backward_sweep = flags & BackwardSweep;
const d Dt = 1 / fs; // Deltat
const d Dt = 1 / _fs; // Deltat
// Estimate N, the number of samples in the sweep part (non-quiescent part):
const us Ns = (us)(Ts * fs);
const us Nq = (us)(Tq * fs);
const us Ns = (us)(Ts * _fs);
const us Nq = (us)(Tq * _fs);
const us N = Ns + Nq;
_signal = vd(N, arma::fill::zeros);
@ -120,7 +133,7 @@ void Sweep::resetImpl() {
d phih = 2 * number_pi * Dt * (fl * Nf + 0.5 * (Nf - 1) * (fu - fl));
us K =
(us)(phih / (2 * number_pi) + Dt * (fu * Nb - (Nb - 1) * (fu - fl)));
(us)(phih / (2 * number_pi) + Dt * (fu * Nb - (Nb - 1) * (fu - fl)));
d eps_num1 = (K - phih / (2 * number_pi)) / Dt;
d eps_num2 = -fu * Nb + (Nb - 1) * (fu - fl);
@ -151,7 +164,7 @@ void Sweep::resetImpl() {
}
} else if (flags & LogSweep) {
DEBUGTRACE_PRINT("Exponential sweep");
DEBUGTRACE_PRINT("Log sweep");
if (forward_sweep || backward_sweep) {
/* Forward or backward sweep */
DEBUGTRACE_PRINT("Forward or backward sweep");
@ -185,9 +198,9 @@ void Sweep::resetImpl() {
const us Nb = N - Nf;
const d k1 = (fu / fl);
const d phif1 =
2 * number_pi * Dt * fl * (k1 - 1) / (d_pow(k1, 1.0 / Nf) - 1);
2 * number_pi * Dt * fl * (k1 - 1) / (d_pow(k1, 1.0 / Nf) - 1);
const us K = (us)(phif1 / (2 * number_pi) +
Dt * fu * (1 / k1 - 1) / (d_pow(1 / k1, 1.0 / Nb) - 1));
Dt * fu * (1 / k1 - 1) / (d_pow(1 / k1, 1.0 / Nb) - 1));
d E;
d k = k1;
@ -196,7 +209,7 @@ void Sweep::resetImpl() {
* continuous */
for (us iter = 0; iter < NITER_NEWTON; iter++) {
E = (k - 1) / (d_pow(k, 1.0 / Nf) - 1) +
(k - 1) / (1 - d_pow(k, -1.0 / Nb)) - K / Dt / fl;
(k - 1) / (1 - d_pow(k, -1.0 / Nb)) - K / Dt / fl;
DEBUGTRACE_PRINT(E);
/* All parts of the derivative of above error E to k */
@ -204,9 +217,9 @@ void Sweep::resetImpl() {
d dEdk2 = (1 / k - 1) / (d_pow(k, -1.0 / Nb) - 1);
d dEdk3 = -1 / (k * (d_pow(k, -1.0 / Nb) - 1));
d dEdk4 = d_pow(k, -1.0 / Nb) * (1 / k - 1) /
(Nb * d_pow(d_pow(k, -1.0 / Nb) - 1, 2));
(Nb * d_pow(d_pow(k, -1.0 / Nb) - 1, 2));
d dEdk5 = -d_pow(k, 1.0 / Nf) * (k - 1) /
(Nf * k * d_pow(d_pow(k, 1.0 / Nf) - 1, 2));
(Nf * k * d_pow(d_pow(k, 1.0 / Nf) - 1, 2));
d dEdk = dEdk1 + dEdk2 + dEdk3 + dEdk4 + dEdk5;
/* Iterate! */

42
src/lasp/dsp/lasp_siggen_impl.h
View File

@ -1,7 +1,19 @@
#pragma once
#include "lasp_siggen.h"
#include "lasp_types.h"
/**
* \addtogroup dsp
* @{
*/
/**
* \addtogroup siggen
* @{
*/
/**
* @brief Generate a random signal (noise)
*/
class Noise : public Siggen {
d level_linear;
virtual vd genSignalUnscaled(const us nframes) override;
@ -19,9 +31,16 @@ class Noise : public Siggen {
};
/**
* @brief Generate a sine wave
*/
class Sine : public Siggen {
d phase = 0;
d omg;
protected:
void resetImpl() override final { phase=0; }
virtual vd genSignalUnscaled(const us nframes) override final;
public:
@ -32,10 +51,14 @@ class Sine : public Siggen {
*/
Sine(const d freq_Hz);
~Sine() = default;
virtual vd genSignalUnscaled(const us nframes) override final;
void setFreq(const d newFreq);
void resetImpl() override final { phase=0; }
void setFreq(const d newFreq) { omg = 2*arma::datum::pi*newFreq; } ;
};
/**
* @brief Base class for all periodic signals (that are exactly periodic based
* on the sampling frequency). Note that the sine wave generator is not exactly
* periodic as the frequency can be any floating point value.
*/
class Periodic: public Siggen {
protected:
vd _signal { 1, arma::fill::zeros};
@ -47,20 +70,23 @@ class Periodic: public Siggen {
};
/**
* @brief Sweep signal
*/
class Sweep : public Periodic {
d fl_, fu_, Ts, Tq;
us index;
us flags;
void resetImpl() override;
public:
static constexpr int ForwardSweep = 1 << 0;
static constexpr int BackwardSweep = 1 << 1;
static constexpr int LinearSweep = 1 << 2;
static constexpr int LogSweep = 1 << 3;
void resetImpl() override;
public:
/**
* Create a sweep signal
*
@ -77,3 +103,5 @@ class Sweep : public Periodic {
~Sweep() = default;
};
/** @} */
/** @} */

2
src/lasp/lasp_cpp.cpp
View File

@ -31,6 +31,7 @@ void init_daqconfiguration(py::module &m);
void init_daq(py::module &m);
void init_streammgr(py::module &m);
void init_datahandler(py::module &m);
void init_siggen(py::module &m);
PYBIND11_MODULE(lasp_cpp, m) {
@ -40,6 +41,7 @@ PYBIND11_MODULE(lasp_cpp, m) {
init_daq(m);
init_streammgr(m);
init_datahandler(m);
init_siggen(m);
// We store the version number of the code via CMake, and create an
// attribute in the C++ code.

9
src/lasp/pybind11/lasp_dsp_pybind.cpp
View File

@ -2,8 +2,6 @@
#include "lasp_avpowerspectra.h"
#include "lasp_biquadbank.h"
#include "lasp_fft.h"
#include "lasp_siggen.h"
#include "lasp_siggen_impl.h"
#include "lasp_slm.h"
#include "lasp_window.h"
#include <iostream>
@ -47,13 +45,6 @@ void init_dsp(py::module &m) {
w.def_static("toTxt", &Window::toText);
/// Siggen
py::class_<Siggen, std::shared_ptr<Siggen>> siggen(m, "Siggen");
siggen.def("setLevel", &Siggen::setLevel,
"Set the level of the signal generator");
py::class_<Sine, std::shared_ptr<Sine>> sw(m, "Sine", siggen);
sw.def(py::init<const d>());
/// SeriesBiquad
py::class_<SeriesBiquad> sbq(m, "SeriesBiquad");

55
src/lasp/pybind11/lasp_siggen.cpp Normal file
View File

@ -0,0 +1,55 @@
#include <carma>
#include "lasp_avpowerspectra.h"
#include "lasp_biquadbank.h"
#include "lasp_fft.h"
#include "lasp_siggen.h"
#include "lasp_siggen_impl.h"
#include "lasp_slm.h"
#include "lasp_window.h"
#include <iostream>
#include <pybind11/pybind11.h>
using std::cerr;
namespace py = pybind11;
/**
* \ingroup pybind
* @{
*
*/
/**
* @brief Initialize Siggen wrappers
*
* @param m The Python module to add classes and methods to
*/
void init_siggen(py::module &m) {
/// Siggen
py::class_<Siggen, std::shared_ptr<Siggen>> siggen(m, "Siggen");
siggen.def("setLevel", &Siggen::setLevel,
"Set the level of the signal generator");
siggen.def("setDCOffset", &Siggen::setDCOffset);
siggen.def("setMute", &Siggen::setMute);
siggen.def("setLevel", &Siggen::setLevel);
siggen.def("setLevel", &Siggen::setLevel, py::arg("newLevel"), py::arg("dB") = true);
siggen.def("genSignal", &Siggen::genSignal);
py::class_<Noise,std::shared_ptr<Noise>> noise(m, "Noise", siggen);
noise.def(py::init<>());
py::class_<Sine, std::shared_ptr<Sine>> sine(m, "Sine", siggen);
sine.def(py::init<const d>());
sine.def("setFreq", &Sine::setFreq);
py::class_<Periodic, std::shared_ptr<Periodic>> periodic(m, "Periodic", siggen);
py::class_<Sweep, std::shared_ptr<Sweep>> sweep(m, "Sweep", periodic);
sweep.def(py::init<const d, const d, const d, const d, const us>());
sweep.def_readonly_static("ForwardSweep", &Sweep::ForwardSweep);
sweep.def_readonly_static("BackwardSweep", &Sweep::BackwardSweep);
sweep.def_readonly_static("LinearSweep", &Sweep::LinearSweep);
sweep.def_readonly_static("LogSweep", &Sweep::LogSweep);
}
/** @} */