From d72a35bfb5a064a9dde2e7c83f4f5bfc0c469dca Mon Sep 17 00:00:00 2001 From: "J.A. de Jong - Redu-Sone B.V., ASCEE V.O.F" Date: Fri, 19 Feb 2021 10:23:51 +0100 Subject: [PATCH] Added logarithmic sweep. For now, chosen to not implement the hyperbolic sweep --- lasp/c/lasp_siggen.c | 600 ++++++++++++++++++++++++++----------------- lasp/c/lasp_siggen.h | 27 +- lasp/wrappers.pyx | 4 + 3 files changed, 392 insertions(+), 239 deletions(-) diff --git a/lasp/c/lasp_siggen.c b/lasp/c/lasp_siggen.c index ba4eaf8..b2d50cf 100644 --- a/lasp/c/lasp_siggen.c +++ b/lasp/c/lasp_siggen.c @@ -5,316 +5,456 @@ // Description: // Signal generator implementation ////////////////////////////////////////////////////////////////////// -#define TRACERPLUS (-5) +/* #define TRACERPLUS (-5) */ #include "lasp_siggen.h" #include "lasp_alloc.h" #include "lasp_assert.h" #include "lasp_mat.h" +/** The fixed number of Newton iterations t.b.d. for tuning the sweep start and + * stop frequency in logarithmic sweeps */ +#define NITER_NEWTON 20 + +/** The number of Bytes of space for the signal-specific data in the Siggen + * structure */ #define PRIVATE_SIZE 64 typedef enum { - SINEWAVE = 0, - NOISE, - SWEEP, - + SINEWAVE = 0, + NOISE, + SWEEP, } SignalType; typedef struct Siggen { - SignalType signaltype; - d fs; // Sampling frequency [Hz] - d level_amp; - char private_data[PRIVATE_SIZE]; + SignalType signaltype; + d fs; // Sampling frequency [Hz] + d level_amp; + char private_data[PRIVATE_SIZE]; } Siggen; typedef struct { - d curtime; - d omg; + d curtime; + d omg; } SinewaveSpecific; typedef struct { - d fl; - d fu; - d Ts; - d phase; - d tau; - bool pos; - us flags; -} SweepSpecific; + us N; + vd data; + us index; +} PeriodicSpecific; typedef struct { - d V1, V2, S; - int phase; - Sosfilterbank* colorfilter; + d V1, V2, S; + int phase; + Sosfilterbank* colorfilter; } NoiseSpecific; static d level_amp(d level_dB){ - return pow(10, level_dB/20); + return pow(10, level_dB/20); } Siggen* Siggen_create(SignalType type, const d fs,const d level_dB) { - fsTRACE(15); + fsTRACE(15); - Siggen* siggen = a_malloc(sizeof(Siggen)); - siggen->signaltype = type; - siggen->fs = fs; - siggen->level_amp = level_amp(level_dB); + Siggen* siggen = a_malloc(sizeof(Siggen)); + siggen->signaltype = type; + siggen->fs = fs; + siggen->level_amp = level_amp(level_dB); - feTRACE(15); - return siggen; + feTRACE(15); + return siggen; } Siggen* Siggen_Sinewave_create(const d fs, const d freq,const d level_dB) { - fsTRACE(15); + fsTRACE(15); - Siggen* sine = Siggen_create(SINEWAVE, fs, level_dB); - dbgassert(sizeof(SinewaveSpecific) <= sizeof(sine->private_data), - "Allocated memory too small"); - SinewaveSpecific* sp = (SinewaveSpecific*) sine->private_data; - sp->curtime = 0; - sp->omg = 2*number_pi*freq; + Siggen* sine = Siggen_create(SINEWAVE, fs, level_dB); + dbgassert(sizeof(SinewaveSpecific) <= sizeof(sine->private_data), + "Allocated memory too small"); + SinewaveSpecific* sp = (SinewaveSpecific*) sine->private_data; + sp->curtime = 0; + sp->omg = 2*number_pi*freq; - feTRACE(15); - return sine; + feTRACE(15); + return sine; } Siggen* Siggen_Noise_create(const d fs, const d level_dB, Sosfilterbank* colorfilter) { - fsTRACE(15); + fsTRACE(15); - Siggen* noise = Siggen_create(NOISE, fs, level_dB); - dbgassert(sizeof(NoiseSpecific) <= sizeof(noise->private_data), - "Allocated memory too small"); - NoiseSpecific* wn = (NoiseSpecific*) noise->private_data; - wn->phase = 0; - wn->V1 = 0; - wn->V2 = 0; - wn->S = 0; - wn->colorfilter = colorfilter; + Siggen* noise = Siggen_create(NOISE, fs, level_dB); + dbgassert(sizeof(NoiseSpecific) <= sizeof(noise->private_data), + "Allocated memory too small"); + NoiseSpecific* wn = (NoiseSpecific*) noise->private_data; + wn->phase = 0; + wn->V1 = 0; + wn->V2 = 0; + wn->S = 0; + wn->colorfilter = colorfilter; - feTRACE(15); - return noise; + feTRACE(15); + return noise; } +Siggen* Siggen_Sweep_create(const d fs,const d fl_,const d fu_, + const d Ts, const us flags, const d level_dB) { + fsTRACE(15); -Siggen* Siggen_Sweep_create(const d fs,const d fl,const d fu, - const d Ts, const us flags, const d level_dB) { - fsTRACE(15); + Siggen* sweep = Siggen_create(SWEEP, fs, level_dB); + dbgassert(sizeof(PeriodicSpecific) <= sizeof(sweep->private_data), + "Allocated memory too small"); - Siggen* sweep = Siggen_create(SWEEP, fs, level_dB); - dbgassert(sizeof(SweepSpecific) <= sizeof(sweep->private_data), - "Allocated memory too small"); - // Set pointer to inplace data storage - SweepSpecific* sp = (SweepSpecific*) sweep->private_data; - if(fl < 0 || fu < 0 || Ts <= 0) { - return NULL; + bool forward_sweep = flags & SWEEP_FLAG_FORWARD; + bool backward_sweep = flags & SWEEP_FLAG_BACKWARD; + + // Set pointer to inplace data storage + dbgassert(!(forward_sweep && backward_sweep), "Both forward and backward flag set"); + + PeriodicSpecific* sp = (PeriodicSpecific*) sweep->private_data; + if(fl_ < 0 || fu_ < 0 || Ts <= 0) { + return NULL; + } + + + const d Dt = 1/fs; // Deltat + + // Estimate N: + const us N = (us) (Ts*fs); + iVARTRACE(15, N); + sp->data = vd_alloc(N); + sp->index = 0; + sp->N = N; + vd* data = &(sp->data); + + // Obtain flags and expand + d phase = 0; + d fl, fu; + + /* Swap fl and fu for a backward sweep */ + if(backward_sweep) { + fu = fl_; + fl = fu_; + } + else { + /* Case of continuous sweep, or forward sweep */ + fl = fl_; + fu = fu_; + } + + /* Linear sweep */ + if(flags & SWEEP_FLAG_LINEAR) { + TRACE(15, "linear sweep"); + + if(forward_sweep || backward_sweep) { + /* Forward or backward sweep */ + TRACE(15, "Forward or backward sweep"); + us K = (us) (Dt*(fl*N+0.5*(N-1)*(fu-fl))); + d eps_num = ((d) K)/Dt - fl*N-0.5*(N-1)*(fu-fl); + d eps = eps_num/(0.5*(N-1)); + iVARTRACE(15, K); + dVARTRACE(15, eps); + + for(us n = 0; n= 0, "BUG"); + + phase += 2*number_pi*Dt*fn; + /* dVARTRACE(17, phase); */ + /* setvecval(data, n, fn); */ + /* setvecval(data, n, phase); */ + } + /* This should be a very small number!! */ + dVARTRACE(15, phase); } - sp->flags = flags; + } + else if(flags & SWEEP_FLAG_EXPONENTIAL) { - sp->fl = fl; - sp->fu = fu; - sp->Ts = Ts; - sp->phase = 0; - sp->pos = flags & SWEEP_FLAG_BACKWARD ? false: true; - if(flags & SWEEP_FLAG_BACKWARD) { - sp->tau = Ts; + TRACE(15, "exponential sweep"); + if(forward_sweep || backward_sweep) { + /* Forward or backward sweep */ + TRACE(15, "Forward or backward sweep"); + d k1 = (fu/fl); + us K = (us) (Dt*fl*(k1-1)/(d_pow(k1,1.0/N)-1)); + d k = k1; + + /* Iterate k to the right solution */ + d E; + for(us iter=0;iter< 10; iter++) { + E = 1 + K/(Dt*fl)*(d_pow(k,1.0/N)-1) - k; + d dEdk = K/(Dt*fl)*d_pow(k,1.0/N)/(N*k)-1; + k -= E/dEdk; + } + + iVARTRACE(15, K); + dVARTRACE(15, k1); + dVARTRACE(15, k); + dVARTRACE(15, E); + + for(us n = 0; ntau = 0; - } - /* sp->pos = false; */ - /* sp->tau = Ts/2; */ + TRACE(15, "Continuous sweep"); - feTRACE(15); - return sweep; + const us Nf = N/2; + 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); + const us K = (us) (phif1/(2*number_pi) + Dt*fu*(1/k1-1)/(d_pow(1/k1,1.0/Nb)-1)); + + d E; + d k = k1; + + /* Newton iterations to converge k to the value such that the sweep is + * continuous */ + for(us iter=0;iter= 0, "BUG"); + + phase += 2*number_pi*Dt*fn; + while(phase > 2*number_pi) phase -= 2*number_pi; + /* dVARTRACE(17, phase); */ + /* setvecval(data, n, fn); */ + /* setvecval(data, n, phase); */ + } + /* This should be a very small number!! */ + dVARTRACE(15, phase); + + } + } + + feTRACE(15); + return sweep; +} + +static void Siggen_periodic_free(PeriodicSpecific* ps) { + assertvalidptr(ps); + fsTRACE(15); + vd_free(&(ps->data)); + feTRACE(15); +} +us Siggen_getN(const Siggen* siggen) { + fsTRACE(15); + assertvalidptr(siggen); + + switch(siggen->signaltype) { + case SINEWAVE: + break; + case NOISE: + break; + case SWEEP: + return ((PeriodicSpecific*) siggen->private_data)->N; + break; + default: + dbgassert(false, "Not implementend signal type"); + + } + + feTRACE(15); + return 0; } void Siggen_free(Siggen* siggen) { - fsTRACE(15); - assertvalidptr(siggen); - NoiseSpecific* sp; + fsTRACE(15); + assertvalidptr(siggen); + NoiseSpecific* sp; - switch(siggen->signaltype) { - case SWEEP: - /* Sweep specific stuff here */ - break; - case SINEWAVE: - /* Sweep specific stuff here */ - break; - case NOISE: - sp = (NoiseSpecific*) siggen->private_data; - if(sp->colorfilter) { - Sosfilterbank_free(sp->colorfilter); - } + switch(siggen->signaltype) { + case SWEEP: + /* Sweep specific stuff here */ + Siggen_periodic_free((PeriodicSpecific*) siggen->private_data); + break; + case SINEWAVE: + /* Sweep specific stuff here */ + break; + case NOISE: + sp = (NoiseSpecific*) siggen->private_data; + if(sp->colorfilter) { + Sosfilterbank_free(sp->colorfilter); + } - } + } - a_free(siggen); - feTRACE(15); + a_free(siggen); + feTRACE(15); } static void Sinewave_genSignal(Siggen* siggen, SinewaveSpecific* sine, vd* samples) { - fsTRACE(15); - assertvalidptr(sine); - d ts = 1/siggen->fs; - d omg = sine->omg; + fsTRACE(14); + assertvalidptr(sine); + d ts = 1/siggen->fs; + d omg = sine->omg; - d curtime = sine->curtime; - for(us i =0; i< samples->n_rows; i++) { - setvecval(samples, i, siggen->level_amp*sin(omg*curtime)); - curtime = curtime + ts; - } - sine->curtime = curtime; - feTRACE(15); + d curtime = sine->curtime; + for(us i =0; i< samples->n_rows; i++) { + setvecval(samples, i, siggen->level_amp*sin(omg*curtime)); + curtime = curtime + ts; + } + sine->curtime = curtime; + feTRACE(14); } -static void Sweep_genSignal(Siggen* siggen, SweepSpecific* sweep, - vd* samples) { - fsTRACE(15); - assertvalidptr(sweep); +static void Periodic_genSignal(Siggen* siggen, PeriodicSpecific* sweep, vd* samples) { + fsTRACE(10); - const d fl = sweep->fl; - const d fu = sweep->fu; - const d deltat = 1/siggen->fs; - const d Ts = sweep->Ts; + for(us i=0; in_rows; i++) { + d* data = getvdval(&(sweep->data), sweep->index); + setvecval(samples, i, siggen->level_amp*(*data)); + sweep->index++; + sweep->index %= sweep->N; + } - const d Thalf = Ts/2; - - dVARTRACE(15, deltat); - - // Load state - d tau = sweep->tau; - bool pos = sweep->pos; - - // Obtain flags and expand - us flags = sweep->flags; - bool forward_sweep = flags & SWEEP_FLAG_FORWARD; - bool backward_sweep = flags & SWEEP_FLAG_BACKWARD; - dbgassert(!(forward_sweep && backward_sweep), "Both forward and backward flag set"); - - d k, Treverse; - if(forward_sweep || backward_sweep) { - k = (fu - fl)/Ts; - Treverse = Ts; - } - else { - k = (fu - fl)/Thalf; - Treverse = Ts/2; - } - - - /* const d k = 0; */ - - d phase = sweep->phase; - d curfreq; - for(us i =0; i< samples->n_rows; i++) { - - curfreq = fl + k*tau; - phase = phase + 2*number_pi*curfreq*deltat; - - // Subtract some to avoid possible overflow. Don't know whether such a - // thing really happens - if(phase > 2*number_pi) - phase = phase - 2*number_pi; - - if(pos) { - tau = tau + deltat; - if(tau >= Treverse) { - if(forward_sweep) { tau = 0; } - else if(backward_sweep) { dbgassert(false, "BUG"); } - else { pos = false; } - } - - } else { - /* dbgassert(false, "cannot get here"); */ - tau = tau - deltat; - if(tau <= 0) { - if(backward_sweep) { tau = Treverse; } - else if(forward_sweep) { dbgassert(false, "BUG"); } - else { pos = true; } - } - } - setvecval(samples, i, siggen->level_amp*d_sin(phase)); - } - // Store state - sweep->phase = phase; - sweep->pos = pos; - sweep->tau = tau; - feTRACE(15); + feTRACE(10); } + static void noise_genSignal(Siggen* siggen, NoiseSpecific* wn, vd* samples) { - fsTRACE(15); - d X; - d S = wn->S; - d V1 = wn->V1; - d V2 = wn->V2; + fsTRACE(15); + d X; + d S = wn->S; + d V1 = wn->V1; + d V2 = wn->V2; - int phase = wn->phase; + int phase = wn->phase; - for(us i =0; i< samples->n_rows; i++) { + for(us i =0; i< samples->n_rows; i++) { - if(wn->phase == 0) { - do { - d U1 = (d)rand() / RAND_MAX; - d U2 = (d)rand() / RAND_MAX; + if(wn->phase == 0) { + do { + d U1 = (d)rand() / RAND_MAX; + d U2 = (d)rand() / RAND_MAX; - V1 = 2 * U1 - 1; - V2 = 2 * U2 - 1; - S = V1 * V1 + V2 * V2; - } while(S >= 1 || S == 0); + V1 = 2 * U1 - 1; + V2 = 2 * U2 - 1; + S = V1 * V1 + V2 * V2; + } while(S >= 1 || S == 0); - X = V1 * sqrt(-2 * d_ln(S) / S); - } else - X = V2 * sqrt(-2 * d_ln(S) / S); + X = V1 * sqrt(-2 * d_ln(S) / S); + } else + X = V2 * sqrt(-2 * d_ln(S) / S); - phase = 1 - phase; + phase = 1 - phase; - setvecval(samples, i, siggen->level_amp*X); - } - if(wn->colorfilter){ - vd filtered = Sosfilterbank_filter(wn->colorfilter, - samples); - dmat_copy(samples, &filtered); - vd_free(&filtered); - } - wn->S = S; - wn->V1 = V1; - wn->V2 = V2; - wn->phase = phase; - feTRACE(15); + setvecval(samples, i, siggen->level_amp*X); + } + if(wn->colorfilter){ + vd filtered = Sosfilterbank_filter(wn->colorfilter, + samples); + dmat_copy(samples, &filtered); + vd_free(&filtered); + } + wn->S = S; + wn->V1 = V1; + wn->V2 = V2; + wn->phase = phase; + feTRACE(15); } void Siggen_genSignal(Siggen* siggen,vd* samples) { - fsTRACE(15); - assertvalidptr(siggen); - assert_vx(samples); + fsTRACE(10); + assertvalidptr(siggen); + assert_vx(samples); - switch(siggen->signaltype) { - case SINEWAVE: - Sinewave_genSignal(siggen, - (SinewaveSpecific*) siggen->private_data, - samples); + switch(siggen->signaltype) { + case SINEWAVE: + Sinewave_genSignal(siggen, + (SinewaveSpecific*) siggen->private_data, + samples); - break; - case NOISE: - noise_genSignal(siggen, - (NoiseSpecific*) siggen->private_data, - samples); - break; - case SWEEP: - Sweep_genSignal(siggen, - (SweepSpecific*) siggen->private_data, - samples); - break; - default: - dbgassert(false, "Not implementend signal type"); + break; + case NOISE: + noise_genSignal(siggen, + (NoiseSpecific*) siggen->private_data, + samples); + break; + case SWEEP: + Periodic_genSignal(siggen, + (PeriodicSpecific*) siggen->private_data, + samples); + break; + default: + dbgassert(false, "Not implementend signal type"); - } + } - feTRACE(15); + feTRACE(10); } diff --git a/lasp/c/lasp_siggen.h b/lasp/c/lasp_siggen.h index e9ed786..2b14b63 100644 --- a/lasp/c/lasp_siggen.h +++ b/lasp/c/lasp_siggen.h @@ -12,6 +12,15 @@ #include "lasp_mat.h" #include "lasp_sosfilterbank.h" +// Define this flag to repeat a forward sweep only, or backward only. If not +// set, we do a continuous sweep +#define SWEEP_FLAG_FORWARD 1 +#define SWEEP_FLAG_BACKWARD 2 + +// Types of sweeps +#define SWEEP_FLAG_LINEAR 4 +#define SWEEP_FLAG_EXPONENTIAL 8 + typedef struct Siggen Siggen; /** @@ -37,15 +46,15 @@ Siggen* Siggen_Sinewave_create(const d fs,const d freq,const d level_dB); */ Siggen* Siggen_Noise_create(const d fs, const d level_dB, Sosfilterbank* colorfilter); -// Define this flag to repeat a forward sweep only, or backward only. If not -// set, we do a continuous sweep -#define SWEEP_FLAG_FORWARD 1 -#define SWEEP_FLAG_BACKWARD 2 -// Types of sweeps -#define SWEEP_FLAG_LINEAR 4 -#define SWEEP_FLAG_EXPONENTIAL 8 -#define SWEEP_FLAG_HYPERBOLIC 16 +/** + * Obtain the repetition period for a periodic excitation. + * @param[in] Siggen* Signal generator handle + * + * @param[out] N The amount of samples in one period, returns 0 if the signal + * does not repeat. + */ +us Siggen_getN(const Siggen*); /** * Create a forward sweep @@ -65,7 +74,7 @@ Siggen* Siggen_Sweep_create(const d fs,const d fl,const d fu, /** * Obtain a new piece of signal * - * @param[in] Siggen* Signal generator private data + * @param[in] Siggen* Signal generator handle * @param[out] samples Samples to fill. Vector should be pre-allocated */ void Siggen_genSignal(Siggen*,vd* samples); diff --git a/lasp/wrappers.pyx b/lasp/wrappers.pyx index b6ee1ea..b7c1163 100644 --- a/lasp/wrappers.pyx +++ b/lasp/wrappers.pyx @@ -635,6 +635,7 @@ cdef extern from "lasp_siggen.h": c_Siggen* Siggen_Sweep_create(d fs, d fl, d fu, d Ts,us sweep_flags, d level_dB) + us Siggen_getN(const c_Siggen*) void Siggen_genSignal(c_Siggen*, vd* samples) nogil void Siggen_free(c_Siggen*) @@ -675,6 +676,9 @@ cdef class Siggen: return output + def getN(self): + return Siggen_getN(self._siggen) + def progress(self): """ TODO: Should be implemented to return the current position in the