From 9caf5fe387bfc07936f022ed29ff629fc4bad683 Mon Sep 17 00:00:00 2001 From: Thijs Hekman Date: Fri, 1 Nov 2024 14:17:31 +0100 Subject: [PATCH] Added amplitude envelope to sweep signals, with associated methods --- cpp_src/dsp/lasp_siggen.h | 2 +- cpp_src/dsp/lasp_siggen_impl.cpp | 131 ++++++++++++++++++++++--------- cpp_src/dsp/lasp_siggen_impl.h | 10 ++- cpp_src/pybind11/lasp_siggen.cpp | 7 ++ 4 files changed, 111 insertions(+), 39 deletions(-) diff --git a/cpp_src/dsp/lasp_siggen.h b/cpp_src/dsp/lasp_siggen.h index 9596efd..fd73bd5 100644 --- a/cpp_src/dsp/lasp_siggen.h +++ b/cpp_src/dsp/lasp_siggen.h @@ -73,7 +73,7 @@ public: * @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 + * @param mute if true */ void setMute(bool mute = true) { _muted = mute; _interruption_frame_count=0; } diff --git a/cpp_src/dsp/lasp_siggen_impl.cpp b/cpp_src/dsp/lasp_siggen_impl.cpp index 7ae0232..f4e7223 100644 --- a/cpp_src/dsp/lasp_siggen_impl.cpp +++ b/cpp_src/dsp/lasp_siggen_impl.cpp @@ -24,60 +24,84 @@ DEBUGTRACE_VARIABLES; Noise::Noise(){DEBUGTRACE_ENTER} -vd Noise::genSignalUnscaled(us nframes) { +vd Noise::genSignalUnscaled(us nframes) +{ return arma::randn(nframes); } void Noise::resetImpl() {} Sine::Sine(const d freq) : omg(2 * arma::datum::pi * freq) { DEBUGTRACE_ENTER; } -vd Sine::genSignalUnscaled(const us nframes) { +vd Sine::genSignalUnscaled(const us nframes) +{ /* DEBUGTRACE_ENTER; */ slock lck(_mtx); const d pi = arma::datum::pi; vd phase_vec = arma::linspace(phase, phase + omg * (nframes - 1) / _fs, nframes); phase += omg * nframes / _fs; - while (phase > 2 * arma::datum::pi) { + while (phase > 2 * arma::datum::pi) + { phase -= 2 * pi; } return arma::sin(phase_vec); } -vd Periodic::genSignalUnscaled(const us nframes) { +vd Periodic::genSignalUnscaled(const us nframes) +{ vd res(nframes); slock lck(_mtx); - if (_signal.size() == 0) { + if (_signal.size() == 0) + { throw rte("No signal defined while calling"); } - for (us i = 0; i < nframes; i++) { - res(i) = _signal[_cur_pos]; + if (_signal.size() != A_.size()) + { + std::cout << "Seq size: " << _signal.size() << ", A size: " << A_.size() << "\n"; + throw rte("Sequence and amplitude envelopes have different lengths"); + } + for (us i = 0; i < nframes; i++) + { + res(i) = A_[_cur_pos] * _signal[_cur_pos]; _cur_pos++; _cur_pos %= _signal.size(); } return res; } +void Periodic::setA(const vd &A) +{ + A_ = A; +} + 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) { + : fl_(fl), fu_(fu), Ts(Ts), Tq(Tq), flags(flags) +{ + if (fl <= 0 || fu < fl || Ts <= 0) + { throw rte("Invalid sweep parameters"); } - if ((flags & ForwardSweep) && (flags & BackwardSweep)) { + if ((flags & ForwardSweep) && (flags & BackwardSweep)) + { throw rte( "Both forward and backward sweep flag set. Please only set either one " "or none for a continuous sweep"); } - if ((flags & LinearSweep) && (flags & LogSweep)) { + if ((flags & LinearSweep) && (flags & LogSweep)) + { throw rte( "Both logsweep and linear sweep flag set. Please only set either one."); } - if (!((flags & LinearSweep) || (flags & LogSweep))) { + if (!((flags & LinearSweep) || (flags & LogSweep))) + { throw rte("Either LinearSweep or LogSweep should be given as flag"); } + + resetImpl(); } -void Sweep::resetImpl() { +void Sweep::resetImpl() +{ DEBUGTRACE_ENTER; slock lck(_mtx); @@ -86,7 +110,7 @@ void Sweep::resetImpl() { 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); @@ -94,14 +118,18 @@ void Sweep::resetImpl() { const us N = Ns + Nq; _signal = vd(N, arma::fill::zeros); + fn_ = vd(N, arma::fill::zeros); index = 0; d fl, fu; /* Swap fl and fu for a backward sweep */ - if (backward_sweep) { + if (backward_sweep) + { fu = fl_; fl = fu_; - } else { + } + else + { /* Case of continuous sweep, or forward sweep */ fl = fl_; fu = fu_; @@ -110,8 +138,10 @@ void Sweep::resetImpl() { d phase = 0; /* Linear sweep */ - if (flags & LinearSweep) { - if (forward_sweep || backward_sweep) { + if (flags & LinearSweep) + { + if (forward_sweep || backward_sweep) + { /* Forward or backward sweep */ /* TRACE(15, "Forward or backward sweep"); */ us K = (us)(Dt * (fl * Ns + 0.5 * (Ns - 1) * (fu - fl))); @@ -120,12 +150,16 @@ void Sweep::resetImpl() { /* iVARTRACE(15, K); */ /* dVARTRACE(15, eps); */ - for (us n = 0; n < Ns; n++) { + for (us n = 0; n < Ns; n++) + { _signal(n) = d_sin(phase); d fn = fl + ((d)n) / Ns * (fu + eps - fl); + fn_(n) = fn; phase += 2 * arma::datum::pi * Dt * fn; } - } else { + } + else + { /* Continous sweep */ /* TRACE(15, "continuous sweep"); */ @@ -150,18 +184,24 @@ void Sweep::resetImpl() { /* dVARTRACE(15, eps); */ d phase = 0; - for (us n = 0; n <= Ns; n++) { + for (us n = 0; n <= Ns; n++) + { /* iVARTRACE(17, n); */ - if (n < N) { + if (n < N) + { _signal[n] = d_sin(phase); } d fn; - if (n <= Nf) { + if (n <= Nf) + { fn = fl + ((d)n) / Nf * (fu - fl); - } else { + } + else + { fn = fu - ((d)n - Nf) / Nb * (fu + eps - fl); } + fn_(n) = fn; /* dbgassert(fn >= 0, "BUG"); */ phase += 2 * number_pi * Dt * fn; @@ -169,9 +209,12 @@ void Sweep::resetImpl() { /* This should be a very small number!! */ /* dVARTRACE(15, phase); */ } - } else if (flags & LogSweep) { + } + else if (flags & LogSweep) + { DEBUGTRACE_PRINT("Log sweep"); - if (forward_sweep || backward_sweep) { + if (forward_sweep || backward_sweep) + { /* Forward or backward sweep */ DEBUGTRACE_PRINT("Forward or backward sweep"); d k1 = (fu / fl); @@ -180,7 +223,8 @@ void Sweep::resetImpl() { /* Iterate k to the right solution */ d E; - for (us iter = 0; iter < 10; iter++) { + for (us iter = 0; iter < 10; iter++) + { E = 1 + K / (Dt * fl) * (d_pow(k, 1.0 / Ns) - 1) - k; d dEdk = K / (Dt * fl) * d_pow(k, 1.0 / Ns) / (Ns * k) - 1; k -= E / dEdk; @@ -191,12 +235,16 @@ void Sweep::resetImpl() { DEBUGTRACE_PRINT(k); DEBUGTRACE_PRINT(E); - for (us n = 0; n < Ns; n++) { + for (us n = 0; n < Ns; n++) + { _signal[n] = d_sin(phase); d fn = fl * d_pow(k, ((d)n) / Ns); + fn_(n) = fn; phase += 2 * number_pi * Dt * fn; } - } else { + } + else + { DEBUGTRACE_PRINT("Continuous sweep"); const us Nf = Ns / 2; @@ -212,7 +260,8 @@ void Sweep::resetImpl() { /* Newton iterations to converge k to the value such that the sweep is * continuous */ - for (us iter = 0; iter < NITER_NEWTON; iter++) { + 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; DEBUGTRACE_PRINT(E); @@ -236,29 +285,37 @@ void Sweep::resetImpl() { DEBUGTRACE_PRINT(k); DEBUGTRACE_PRINT(E); - for (us n = 0; n <= Ns; n++) { + for (us n = 0; n <= Ns; n++) + { /* iVARTRACE(17, n); */ - if (n < Ns) { + if (n < Ns) + { _signal[n] = d_sin(phase); } d fn; - if (n <= Nf) { + if (n <= Nf) + { fn = fl * d_pow(k, ((d)n) / Nf); - } else { + } + else + { fn = fl * k * d_pow(1 / k, ((d)n - Nf) / Nb); } + fn_(n) = fn; /* dbgassert(fn >= 0, "BUG"); */ phase += 2 * number_pi * Dt * fn; - while (phase > 2 * number_pi) phase -= 2 * number_pi; + while (phase > 2 * number_pi) + phase -= 2 * number_pi; /* dVARTRACE(17, phase); */ } /* This should be a very small number!! */ DEBUGTRACE_PRINT(phase); } - } // End of log sweep - else { + } // End of log sweep + else + { // Either log or linear sweep had to be given as flags. assert(false); } diff --git a/cpp_src/dsp/lasp_siggen_impl.h b/cpp_src/dsp/lasp_siggen_impl.h index 9069b1e..e5847cc 100644 --- a/cpp_src/dsp/lasp_siggen_impl.h +++ b/cpp_src/dsp/lasp_siggen_impl.h @@ -58,6 +58,7 @@ class Sine : public Siggen { * periodic as the frequency can be any floating point value. */ class Periodic: public Siggen { + protected: vd _signal { 1, arma::fill::zeros}; us _cur_pos = 0; @@ -68,8 +69,12 @@ class Periodic: public Siggen { * @return As stated above */ vd getSequence() const { return _signal; } + vd A_ { 1, arma::fill::ones}; + + void setA(const vd& A); virtual vd genSignalUnscaled(const us nframes) override final; + ~Periodic() = default; }; @@ -81,7 +86,8 @@ class Sweep : public Periodic { d fl_, fu_, Ts, Tq; us index; us flags; - + vd fn_ { 1, arma::fill::zeros}; + void resetImpl() override; public: @@ -90,6 +96,8 @@ class Sweep : public Periodic { static constexpr int LinearSweep = 1 << 2; static constexpr int LogSweep = 1 << 3; + vd getfn() const { return fn_; } + /** * Create a sweep signal * diff --git a/cpp_src/pybind11/lasp_siggen.cpp b/cpp_src/pybind11/lasp_siggen.cpp index 246133b..0161069 100644 --- a/cpp_src/pybind11/lasp_siggen.cpp +++ b/cpp_src/pybind11/lasp_siggen.cpp @@ -43,11 +43,18 @@ void init_siggen(py::module &m) { py::class_> periodic(m, "Periodic", siggen); + periodic.def("setA", + [](Periodic &p, const dpyarray A) { + p.setA(NpyToCol(A)); + }); + periodic.def("getSequence", [](const Sweep &s) { return ColToNpy(s.getSequence()); }); py::class_> sweep(m, "Sweep", periodic); sweep.def(py::init()); + sweep.def("getfn", + [](const Sweep &s) { return ColToNpy(s.getfn()); }); sweep.def_readonly_static("ForwardSweep", &Sweep::ForwardSweep); sweep.def_readonly_static("BackwardSweep", &Sweep::BackwardSweep); sweep.def_readonly_static("LinearSweep", &Sweep::LinearSweep);