API — ICoreBlocks/ICoreMath/SignalProcessing
The public contract of 24 header(s) under src/ICoreBlocks/ICoreMath/SignalProcessing — 24 class/struct definition(s), 18 declaration(s). Each section shows the header's banner and its public (and protected-virtual) surface exactly as the file writes it.
ICoreAnalogPrototypes.h#
src/ICoreBlocks/ICoreMath/SignalProcessing/ICoreAnalogPrototypes.h
ICoreAnalogPrototypes#
ICoreAnalogPrototypes.h:33 · class · nested Prototype · 2 declaration(s)
The five analog lowpass prototypes every IIR design starts from (toolbox row T2.5): each answers a zero set, a pole set and a gain for a filter whose passband edge is 1 rad/s, which the band transf...
class ICoreAnalogPrototypes {
public:
// What every prototype answers: a zero set (empty for three of the five),
// a pole set, and the gain that makes the passband what the approximation
// promises.
struct Prototype {
std::vector<std::complex<double>> zeros;
std::vector<std::complex<double>> poles;
double gain = 1.0;
};
// MATLAB's buttap(n): n poles evenly spaced on the left half of the unit
// circle, no zeros, gain 1. The odd order's real pole at -1 comes LAST,
// which is MATLAB's order and is reproduced.
static bool butterworth(const size_t& order, Prototype& out, std::string* whyNot = nullptr);
// MATLAB's cheb1ap(n, rp): the Butterworth circle squashed onto an ellipse
// by the passband ripple `rp` in dB. An EVEN order does not reach unit
// gain at DC -- its gain is patched down by sqrt(1 + epsilon^2), which is
// MATLAB's rule and the reason an even Chebyshev I starts a ripple low.
static bool chebyshevI(const size_t& order, const double& passbandRipple,
Prototype& out, std::string* whyNot = nullptr);
// MATLAB's cheb2ap(n, rs): the RECIPROCAL of a Chebyshev I geometry, so
// the ripple lands in the stopband and the filter gains zeros on the
// imaginary axis. `rs` is the stopband attenuation in dB.
static bool chebyshevII(const size_t& order, const double& stopbandAttenuation,
Prototype& out, std::string* whyNot = nullptr);
// MATLAB's ellipap(n, rp, rs), by way of its `ellipap2`: the poles and
// zeros are values of the Jacobi elliptic functions at the odd quarter
// periods, and the modulus that ties the two ripples to the order is the
// solution of the degree equation.
//
// `rp` must be above zero and below `rs` -- an elliptic filter with no
// passband ripple is a Chebyshev II, and one whose passband ripple exceeds
// its stopband attenuation is not a filter at all. Both are refused with
// that reason, as MATLAB refuses them.
static bool elliptic(const size_t& order, const double& passbandRipple,
const double& stopbandAttenuation, Prototype& out,
std::string* whyNot = nullptr);
// MATLAB's besselap(n): the poles are the roots of the reverse Bessel
// polynomial, whose coefficients are exact integers,
//
// theta_n(s) = sum_k [ (2n-k)! / (2^(n-k) k! (n-k)!) ] s^k
//
// scaled by the n-th root of that polynomial's constant term. **That
// scaling is the whole of MATLAB's convention and it is not the -3 dB
// one**: it makes prod(-p) exactly 1, so the gain is 1 and the DC gain is
// 1 with a monic denominator, and the filter's -3 dB point is therefore
// NOT at 1 rad/s (measured: 1.756 rad/s at order 3, 2.114 at order 4).
// Normalising to the -3 dB frequency instead -- the obvious reading of
// "prototype" -- answers poles 2.4 times too large at order 3.
//
// MATLAB refuses an order above 25 (its own table stops there); this
// computes the roots and so has no table, but the refusal is kept because
// the polynomial's coefficients overflow a double soon after.
static bool bessel(const size_t& order, Prototype& out, std::string* whyNot = nullptr);
};
};ICoreAnalyticSignal.h#
src/ICoreBlocks/ICoreMath/SignalProcessing/ICoreAnalyticSignal.h
ICoreAnalyticSignal#
ICoreAnalyticSignal.h:43 · class · nested Envelope · 0 declaration(s)
The analytic signal and the envelopes read off it (toolbox row T2.19): hilbert and envelope.
class ICoreAnalyticSignal {
public:
// MATLAB's hilbert(x[, n]). `length` of 0 asks for the signal's own
// length; any other value is MATLAB's `n`, and the answer has that many
// samples. The orientation of the input is kept.
static bool analytic(const ICoreMatrix& signal, const size_t& length,
ICoreComplexMatrix& out, std::string* whyNot = nullptr);
// Which envelope is asked for. `Analytic` with a window length is
// MATLAB's FIR form; without one it is the transform above.
enum class Method { Analytic, Rms, Peaks };
// The two envelopes, which travel together because a lower envelope on its
// own is a curve nobody can place.
struct Envelope {
ICoreMatrix upper;
ICoreMatrix lower;
};
// MATLAB's [yupper, ylower] = envelope(x[, n, method]). `window` is
// MATLAB's n: the FIR order for `Analytic`, the sliding length for `Rms`
// and the minimum peak separation for `Peaks`; 0 asks for the plain
// analytic form and is refused by the other two, which have no default.
static bool envelope(const ICoreMatrix& signal, const Method& method,
const size_t& window, Envelope& out, std::string* whyNot = nullptr);
};
};ICoreCorrelation.h#
src/ICoreBlocks/ICoreMath/SignalProcessing/ICoreCorrelation.h
ICoreCorrelation#
ICoreCorrelation.h:31 · class · 0 declaration(s)
The Signal Processing Toolbox correlation family that base MATLAB's xcorr leaves out (toolbox row T2.13): the covariance spelling, the 2-D correlation, circular convolution, and the delay estimat...
class ICoreCorrelation {
public:
// MATLAB's xcov: ICoreMatrix::xcorr of the MEAN-REMOVED signals. The
// arguments mean exactly what xcorr's do -- `y` NULL asks for the
// autocovariance, `maxlag` below zero asks for MATLAB's default
// max(numel(x), numel(y)) - 1, `scale` is "none" (the default), "biased",
// "unbiased" or "coeff" -- and the orientation follows `x`.
//
// The mean removal happens FIRST, so "coeff" normalizes by the
// mean-removed energies and xcov(x, x, 'coeff') is 1 at lag 0 for any
// signal with a non-zero variance, which is the property that makes the
// covariance spelling worth having.
static ICoreMatrix xcov(const ICoreMatrix& x, const ICoreMatrix* y,
const int& maxlag, const std::string& scale,
ICoreMatrix* lags, std::string* whyNot = nullptr);
// MATLAB's xcorr2: the 2-D cross-correlation of two matrices,
//
// C(p, q) = sum_{i, j} A(i, j) * B(i - p, j - q)
//
// over every 2-D lag, so the answer is (ma + mb - 1) x (na + nb - 1) with
// lag (0, 0) at its centre entry (mb, nb) counting from one. `b` NULL asks
// for the 2-D AUTOcorrelation, which is MATLAB's one-argument form.
static ICoreMatrix xcorr2(const ICoreMatrix& a, const ICoreMatrix* b,
std::string* whyNot = nullptr);
// MATLAB's cconv(a, b, n): the length-`n` CIRCULAR convolution. Both
// sequences are wrapped modulo `n` first (MATLAB's `datawrap`, so a
// sequence longer than `n` folds onto itself rather than being cut), and
// the result is the circular convolution of the two wrapped sequences.
// `n` = 0 asks for MATLAB's default, numel(a) + numel(b) - 1, at which
// length the circular convolution IS the linear one.
//
// This is the direct O(n^2) sum where MATLAB's is an FFT, an elementwise
// product and an inverse FFT, so the two agree to rounding rather than
// exactly -- and, as with xcorr, it is MATLAB that carries the error.
static ICoreMatrix cconv(const ICoreMatrix& a, const ICoreMatrix& b,
const size_t& n, std::string* whyNot = nullptr);
// MATLAB's finddelay(x, y[, maxlag]): the estimated delay of `y` relative
// to `x`, in samples, as a whole number in [-maxlag, maxlag]. `maxlag`
// below zero asks for MATLAB's default, max(numel(x), numel(y)) - 1.
//
// Positive means `y` is LATER than `x`. The estimate is the lag of the
// largest |xcorr| normalized by sqrt(sum(x.^2) * sum(y.^2)); a tie between
// the best positive and the best negative lag goes to the smaller |lag|,
// and to the POSITIVE one when the two are the same size; and a peak below
// 1e-8 answers 0, because a correlation that small does not name a delay.
static double findDelay(const ICoreMatrix& x, const ICoreMatrix& y,
const int& maxlag, std::string* whyNot = nullptr);
// MATLAB's [xa, ya, d] = alignsignals(x, y[, maxlag[, 'truncate']]):
// `d` is findDelay's answer, and the signal that starts EARLY is delayed
// by |d| leading zeros so the two line up. With `truncate` the padded
// signal keeps its original length (its tail falls off the end) instead of
// growing by |d|; a delay at least as large as the signal then leaves all
// zeros, which is MATLAB's documented behaviour and its warning.
//
// The orientation of each output follows its own input.
static bool alignSignals(const ICoreMatrix& x, const ICoreMatrix& y,
const int& maxlag, const bool& truncate,
ICoreMatrix& xa, ICoreMatrix& ya, double& delay,
std::string* whyNot = nullptr);
};
};ICoreDiscreteTransforms.h#
src/ICoreBlocks/ICoreMath/SignalProcessing/ICoreDiscreteTransforms.h
ICoreDiscreteTransforms#
ICoreDiscreteTransforms.h:31 · class · 0 declaration(s)
The discrete transforms beside the Fourier one (toolbox row T2.20): the cosine family, and the two ways to evaluate a spectrum somewhere other than on the FFT's own grid.
class ICoreDiscreteTransforms {
public:
// MATLAB's dct / idct. `type` is 1..4 (2 is MATLAB's default), `length` is
// the transform length -- 0 asks for the signal's own, and any other value
// truncates or zero-pads first, as MATLAB's `dct(x, n)` does. `inverse`
// asks for `idct`, which is the same matrix transposed.
//
// A MATRIX is transformed COLUMN BY COLUMN (MATLAB's rule); a vector keeps
// its orientation.
static ICoreMatrix cosineTransform(const ICoreMatrix& x, const size_t& length,
const int& type, const bool& inverse,
std::string* whyNot = nullptr);
// MATLAB's czt(x, m, w, a): the chirp z-transform, which evaluates
//
// X(j) = sum_n x(n) * (a * w^-j)^-n
//
// at `m` points along the spiral z_j = a * w^(-j). The defaults are
// MATLAB's -- m = numel(x), w = exp(-2i*pi/m), a = 1 -- at which the
// answer IS the DFT.
//
// This is the direct O(m*n) sum where MATLAB's is Bluestein's algorithm
// (three FFTs), so the two agree to rounding rather than exactly, and it is
// MATLAB carrying the error: its czt of [1 2 3 4] answers a 4.6e-17
// imaginary part where the sum is real.
static ICoreComplexMatrix chirpZ(const ICoreMatrix& x, const size_t& points,
const std::complex<double>& ratio,
const std::complex<double>& start,
std::string* whyNot = nullptr);
// MATLAB's goertzel(x, k): the DFT at the frequency indices `k`, counted
// from ONE as MATLAB counts them, so k = 1 is the DC bin. `indices` empty
// asks for every bin, which is `fft`.
//
// The indices need not be whole numbers: MATLAB evaluates
// X(k) = sum_n x(n) * exp(-2i*pi*(k-1)*(n-1)/N) for any real k, and a
// fractional one lands between bins rather than being refused.
static ICoreComplexMatrix goertzel(const ICoreMatrix& x, const ICoreMatrix* indices,
std::string* whyNot = nullptr);
};
};ICoreEllipticFunctions.h#
src/ICoreBlocks/ICoreMath/SignalProcessing/ICoreEllipticFunctions.h
ICoreEllipticFunctions#
ICoreEllipticFunctions.h:31 · class · 0 declaration(s)
The complete elliptic integral and the three Jacobi elliptic functions an elliptic filter is built out of (toolbox rows T2.5 and T2.3).
class ICoreEllipticFunctions {
public:
// The descending Landen sequence of moduli, MATLAB's `landen(k, tol)`:
// each entry is the next modulus down, stopping once one falls at or below
// `tolerance` (machine epsilon by default). k = 0 or 1 answers itself.
static std::vector<double> landen(const double& modulus,
const double& tolerance = 2.220446049250313e-16);
// The complete elliptic integral K(k) and its complement K'(k) = K(k'),
// MATLAB's `ellipk`. Either may be infinite -- K(1) and K'(0) are -- and
// that is an answer rather than a failure.
static void completeIntegral(const double& modulus, double& integral, double& complement,
const double& tolerance = 2.220446049250313e-16);
// The Jacobi elliptic functions cd() and sn() in MATLAB's normalized
// spelling, where the argument is measured in QUARTER PERIODS: `cde(u, k)`
// is cd(u*K(k), k) and `sne(u, k)` is sn(u*K(k), k). Both take a complex
// argument, which is what an elliptic filter's poles need.
static std::complex<double> cde(const std::complex<double>& u, const double& modulus,
const double& tolerance = 2.220446049250313e-16);
static std::complex<double> sne(const std::complex<double>& u, const double& modulus,
const double& tolerance = 2.220446049250313e-16);
// Their inverses, again in quarter periods: `acde(w, k)` solves
// cde(u, k) = w and `asne(w, k)` solves sne(u, k) = w, both by the
// ASCENDING Landen transformation. The answer is brought into
// -2 < Re(u) <= 2 and |Im(u)| <= K'/K, as MATLAB's is.
static std::complex<double> acde(const std::complex<double>& w, const double& modulus,
const double& tolerance = 2.220446049250313e-16);
static std::complex<double> asne(const std::complex<double>& w, const double& modulus,
const double& tolerance = 2.220446049250313e-16);
// The DEGREE equation, MATLAB's `ellipdeg(N, k1)`: the selectivity modulus
// k of an order-N elliptic filter whose ripple modulus is k1 -- the
// solution of N * K'(k)/K(k) = K'(k1)/K(k1). Below k1 = 1e-6 it is solved
// through the nome instead, because the Landen sequence loses the
// difference between k1 and 0 there.
static double degree(const size_t& order, const double& rippleModulus,
const double& tolerance = 2.220446049250313e-16);
// The degree equation the other way, MATLAB's `ellipdeg2(1/N, k)`: the
// ripple modulus whose order-N filter has selectivity modulus k. `order`
// is the RECIPROCAL exponent, so pass 1/N to answer k1 from k -- which is
// the direction `ellipord` needs.
static double degreeInverse(const double& exponent, const double& modulus,
const double& tolerance = 2.220446049250313e-16);
};
};ICoreFFT.h#
src/ICoreBlocks/ICoreMath/SignalProcessing/ICoreFFT.h
ICoreFFT#
ICoreFFT.h:8 · class · 4 declaration(s)
Discrete Fourier transform of a logged/simulated signal, backed by Eigen::FFT (kissfft backend, bundled -- no external dependency).
class ICoreFFT {
public:
// signal: a row or column vector, N samples.
// Returns the complex spectrum as an N x 2 matrix: column 0 = real part,
// column 1 = imaginary part.
static ICoreMatrix forward(const ICoreMatrix& signal);
// spectrum: an N x 2 matrix [real, imag], as returned by forward().
// Returns the reconstructed N x 1 time-domain signal.
static ICoreMatrix inverse(const ICoreMatrix& spectrum);
// Magnitude / phase (radians) of forward(signal), each as an N x 1 vector.
static ICoreMatrix magnitude(const ICoreMatrix& signal);
static ICoreMatrix phase(const ICoreMatrix& signal);
};
};ICoreFilterAnalysis.h#
src/ICoreBlocks/ICoreMath/SignalProcessing/ICoreFilterAnalysis.h
ICoreFilterAnalysis#
ICoreFilterAnalysis.h:34 · class · 0 declaration(s)
What a designed filter is LOOKED AT with in the time domain and in phase (toolbox row T2.7): its impulse and step responses, its group delay and its unwrapped phase.
class ICoreFilterAnalysis {
public:
// MATLAB's impzlength(b, a): how long the response has to be for the
// filter to have finished. `tolerance` is MATLAB's 5e-5.
static bool responseLength(const ICoreMatrix& numerator, const ICoreMatrix& denominator,
const double& tolerance, size_t& out,
std::string* whyNot = nullptr);
// MATLAB's [h, t] = impz(b, a[, n][, fs]) and [s, t] = stepz(...): the same
// function with a different input. `points` of 0 asks for responseLength();
// a non-null `indices` is MATLAB's vector form, which answers the response
// AT those sample indices and may reach back before zero.
static bool timeResponse(const ICoreMatrix& numerator, const ICoreMatrix& denominator,
const ICoreMatrix* indices, const size_t& points,
const double& sampleRate, const bool& stepInput,
ICoreMatrix& values, ICoreMatrix& times,
std::string* whyNot = nullptr);
// MATLAB's [gd, w] = grpdelay(b, a[, n][, 'whole'][, fs]): the negative
// derivative of the phase, in SAMPLES. `frequencies` is the caller's own
// grid (in rad/sample when `sampleRate` is 0, in Hz otherwise); otherwise
// `points` picks the uniform grid, 0 asking for MATLAB's default of 512.
static bool groupDelay(const ICoreMatrix& numerator, const ICoreMatrix& denominator,
const ICoreMatrix* frequencies, const size_t& points,
const bool& whole, const double& sampleRate, ICoreMatrix& delay,
ICoreMatrix& grid, std::string* whyNot = nullptr);
// MATLAB's [phi, w] = phasez(b, a[, n][, 'whole'][, fs]): the UNWRAPPED
// phase, which is the one an analysis reads -- `angle(freqz(...))` is the
// same curve folded into (-pi, pi] and is not what this answers.
static bool unwrappedPhase(const ICoreMatrix& numerator, const ICoreMatrix& denominator,
const size_t& points, const bool& whole, const double& sampleRate,
ICoreMatrix& phase, ICoreMatrix& grid,
std::string* whyNot = nullptr);
};
};ICoreFilterConversions.h#
src/ICoreBlocks/ICoreMath/SignalProcessing/ICoreFilterConversions.h
ICoreFilterConversions#
ICoreFilterConversions.h:43 · class · nested StateSpace, ZeroPoleGain · 0 declaration(s)
The four ways a filter can be written down, and the conversions between them (toolbox row T2.8): coefficients, zeros-poles-gain, state space, and a cascade of second-order sections.
class ICoreFilterConversions {
public:
// A realization. `a` is nx by nx, `b` nx by 1, `c` ny by nx, `d` ny by 1.
struct StateSpace {
ICoreMatrix a;
ICoreMatrix b;
ICoreMatrix c;
ICoreMatrix d;
};
// Zeros, poles and gain. `zeros` is nz by ny (Inf-padded when the outputs
// have different zero counts, which is MATLAB's own shape), `poles` a
// column, `gain` ny by 1.
struct ZeroPoleGain {
ICoreComplexMatrix zeros;
ICoreComplexMatrix poles;
ICoreMatrix gain;
};
// zp2sos's section ORDER: 'up' puts the pole pair furthest from the unit
// circle first, 'down' reverses the cascade. MATLAB's default is Up.
enum class SectionOrder { Up, Down };
// MATLAB's [A, B, C, D] = tf2ss(num, den): the CONTROLLER canonical form.
// `num` may have one row per output; `den` is a row.
static bool tfToStateSpace(const ICoreMatrix& numerator, const ICoreMatrix& denominator,
StateSpace& out, std::string* whyNot = nullptr);
// MATLAB's [z, p, k] = tf2zp(num, den). Leading zero columns of the
// numerator are stripped first -- they are not zeros of the filter, they
// are absent leading terms -- and `k` is the first surviving numerator
// coefficient over the denominator's leading one.
static bool tfToZeroPoleGain(const ICoreMatrix& numerator, const ICoreMatrix& denominator,
ZeroPoleGain& out, std::string* whyNot = nullptr);
// MATLAB's [num, den] = zp2tf(z, p, k). Both polynomials are the REAL part
// of the product over the roots, which is what makes a conjugate-paired
// set answer real coefficients rather than a complex vector with zero
// imaginary parts.
static bool zeroPoleGainToTf(const ICoreComplexMatrix& zeros, const ICoreComplexMatrix& poles,
const ICoreMatrix& gain, ICoreMatrix& numerator,
ICoreMatrix& denominator, std::string* whyNot = nullptr);
// MATLAB's [z, p, k] = ss2zp(A, B, C, D, iu): the response from the iu'th
// input (1-based, as MATLAB counts). Poles are eig(A) for every output;
// zeros and gain are per output row.
static bool stateSpaceToZeroPoleGain(const ICoreMatrix& a, const ICoreMatrix& b,
const ICoreMatrix& c, const ICoreMatrix& d,
const size_t& input, ZeroPoleGain& out,
std::string* whyNot = nullptr);
// MATLAB's [num, den] = ss2tf(A, B, C, D, iu), which is ss2zp followed by
// poly() on each side -- the route matters, because it is what makes
// ss2tf's numerator carry the zeros' rounding rather than the state
// matrices'.
static bool stateSpaceToTf(const ICoreMatrix& a, const ICoreMatrix& b, const ICoreMatrix& c,
const ICoreMatrix& d, const size_t& input, ICoreMatrix& numerator,
ICoreMatrix& denominator, std::string* whyNot = nullptr);
// MATLAB's [A, B, C, D] = zp2ss(z, p, k): the BALANCED cascade of
// second-order sections, not the controller form -- pairs of poles are
// realized together and each section is balanced by the geometric mean of
// its pole magnitudes, which is why this and tf2ss answer different
// matrices for the same filter and both are right.
static bool zeroPoleGainToStateSpace(const ICoreComplexMatrix& zeros,
const ICoreComplexMatrix& poles, const ICoreMatrix& gain,
StateSpace& out, std::string* whyNot = nullptr);
// MATLAB's [sos, g] = zp2sos(z, p, k, order). The gain is separated only
// when the caller asks for it; the single-output form folds it into the
// FIRST section's numerator, which this reports through `gain` and leaves
// to the caller to fold (the evaluator does, so both console forms match
// MATLAB's two).
static bool zeroPoleGainToSections(const ICoreComplexMatrix& zeros,
const ICoreComplexMatrix& poles, const ICoreMatrix& gain,
const SectionOrder& order, ICoreMatrix& sections,
double& sectionGain, std::string* whyNot = nullptr);
// MATLAB's [sos, g] = tf2sos(b, a) and ss2sos(A, B, C, D, iu): the same
// pairing, reached from the other two spellings. tf2sos equalizes the
// coefficient lengths first (MATLAB's `eqtflength`), which is what makes
// tf2sos([1 2], [1 3 2]) a two-pole cascade rather than an error.
static bool tfToSections(const ICoreMatrix& numerator, const ICoreMatrix& denominator,
const SectionOrder& order, ICoreMatrix& sections,
double& sectionGain, std::string* whyNot = nullptr);
static bool stateSpaceToSections(const ICoreMatrix& a, const ICoreMatrix& b,
const ICoreMatrix& c, const ICoreMatrix& d,
const size_t& input, const SectionOrder& order,
ICoreMatrix& sections, double& sectionGain,
std::string* whyNot = nullptr);
// MATLAB's [b, a] = sos2tf(sos, g) and [z, p, k] = sos2zp(sos, g): the
// cascade multiplied back out. A section whose third and sixth entries are
// both zero is a FIRST-order section and its trailing zeros are dropped
// before its roots are taken -- otherwise every such section would
// contribute a spurious zero at the origin.
static bool sectionsToTf(const ICoreMatrix& sections, const ICoreMatrix& gain,
ICoreMatrix& numerator, ICoreMatrix& denominator,
std::string* whyNot = nullptr);
static bool sectionsToZeroPoleGain(const ICoreMatrix& sections, const ICoreMatrix& gain,
ZeroPoleGain& out, std::string* whyNot = nullptr);
// MATLAB's cplxpair, exposed because zp2sos and zp2ss both need it and
// because it is a name of its own (core board C10.10): conjugate pairs
// adjacent with the negative imaginary part first, ordered by ascending
// real part, and the real entries sorted at the END. A set that does not
// pair within `tolerance` is refused rather than silently reordered --
// which is the behaviour zp2ss depends on when it retries with a looser
// tolerance.
static bool complexPair(std::vector<std::complex<double>>& values, const double& tolerance,
std::string* whyNot = nullptr);
};
};ICoreFilterDesign.h#
src/ICoreBlocks/ICoreMath/SignalProcessing/ICoreFilterDesign.h
ICoreFilterDesign#
ICoreFilterDesign.h:9 · class · 2 declaration(s)
Digital Butterworth filter design via the analog prototype + bilinear (Tustin) transform, with pre-warped band edges.
class ICoreFilterDesign {
public:
// order: filter order (>= 1). cutoffHz: -3 dB cutoff frequency, must be
// in (0, sampleRateHz/2). sampleRateHz: discrete sampling rate.
// Returns a discrete (Ts = 1/sampleRateHz) transfer function normalized
// to unity gain in the passband (DC for lowpass, Nyquist for highpass).
static ICoreTransferFunction butterworthLowpass(const size_t& order, const double& cutoffHz, const double& sampleRateHz);
static ICoreTransferFunction butterworthHighpass(const size_t& order, const double& cutoffHz, const double& sampleRateHz);
// Band-stop (notch) of the given prototype order, built with the
// lowpass-to-bandstop transform s -> B*s / (s^2 + w0^2). It nulls
// centerHz and passes everything either side of it.
//
// bandwidthHz is the -3 dB width of the stop band, and centerHz +/-
// bandwidthHz/2 must both fall inside (0, sampleRateHz/2). Each prototype
// pole becomes two, so the result has 2*order poles and 2*order zeros --
// order 1 is the classic two-pole/two-zero notch, and higher orders square
// up the shoulders without moving the null.
//
// The zeros land exactly on the unit circle at the centre frequency, so the
// attenuation there is total rather than merely deep (measured |H| at the
// centre is ~1e-13 for order 1, ~1e-10 for order 3, the residual being
// coefficient round-off from the repeated roots). Gain is normalized to
// unity at DC, which for a band-stop is in the passband.
//
// Note on the band edges: this filter shape is geometrically symmetric about
// the centre, which is two degrees of freedom -- not enough to place the null
// at centerHz AND both -3 dB points at centerHz +/- bandwidthHz/2. The null
// is kept exact, so the achieved -3 dB edges come out geometrically
// symmetric (sqrt(f_lo*f_hi) = centerHz) rather than arithmetically. Their
// separation still matches bandwidthHz -- within 0.01% for a narrow notch,
// and a few tenths of a percent for one spanning a large fraction of the
// band -- so the width you ask for is the width you get; only its placement
// shifts slightly downward relative to an arithmetic split.
static ICoreTransferFunction butterworthNotch(const size_t& order, const double& centerHz,
const double& bandwidthHz, const double& sampleRateHz);
};
};ICoreFilterTransforms.h#
src/ICoreBlocks/ICoreMath/SignalProcessing/ICoreFilterTransforms.h
ICoreFilterTransforms#
ICoreFilterTransforms.h:31 · class · 0 declaration(s)
The frequency and s-to-z transformations a filter design is assembled from (toolbox row T2.9): move an analog prototype to the band you want, then carry it to the discrete plane.
class ICoreFilterTransforms {
public:
// Which band the analog prototype is being moved to.
enum class Band {
Lowpass, // lp2lp: s -> s / Wo
Highpass, // lp2hp: s -> Wo / s
Bandpass, // lp2bp: s -> (s^2 + Wo^2) / (Bw * s)
Bandstop, // lp2bs: s -> (Bw * s) / (s^2 + Wo^2)
};
// MATLAB's lp2lp / lp2hp / lp2bp / lp2bs on a coefficient PAIR: the
// prototype num/den in descending powers of s, and the answer in the same
// form with the denominator made monic -- which is what MATLAB's
// `poly(at)` route produces and is therefore part of the contract, not
// tidying.
//
// `bandwidth` is read only for the two band forms. A transform that leaves
// no leading denominator coefficient to divide by -- lp2hp of a prototype
// with a pole at the origin -- is refused with that reason rather than
// answered as an infinity.
static bool bandTransform(const ICoreMatrix& numerator, const ICoreMatrix& denominator,
const Band& band, const double& cutoff, const double& bandwidth,
ICoreMatrix& outNumerator, ICoreMatrix& outDenominator,
std::string* whyNot = nullptr);
// MATLAB's [bd, ad] = bilinear(b, a, fs[, fp]): the analog pair carried to
// the discrete plane by s = fs2*(z-1)/(z+1), with fs2 = 2*fs, and the
// answer made monic.
//
// `prewarp` is MATLAB's optional fp in Hz: it replaces fs2 with
// 2*pi*fp / tan(pi*fp/fs), which is what makes the analog and digital
// responses agree AT fp rather than only near zero. Pass a null pointer
// for the unwarped transform.
//
// The numerator may not be of higher degree than the denominator, which is
// MATLAB's own restriction (an improper analog filter has no bilinear
// image of the same order).
static bool bilinear(const ICoreMatrix& numerator, const ICoreMatrix& denominator,
const double& sampleRate, const double* prewarp,
ICoreMatrix& outNumerator, ICoreMatrix& outDenominator,
std::string* whyNot = nullptr);
// MATLAB's [zd, pd, kd] = bilinear(z, p, k, fs[, fp]): the same transform
// on a zero-pole-gain triple, which maps each zero and pole by
// (1 + s/fs2) / (1 - s/fs2) and scales the gain by
// prod(fs2 - z) / prod(fs2 - p).
//
// The two shapes MATLAB's own answer has, and both matter: the zeros come
// back PADDED with as many -1 entries as the filter has excess poles (the
// bilinear map sends infinity to z = -1), and an infinite zero on the way
// in is dropped before the count is taken.
static bool bilinearZeroPoleGain(const ICoreComplexMatrix& zeros,
const ICoreComplexMatrix& poles, const double& gain,
const double& sampleRate, const double* prewarp,
ICoreComplexMatrix& outZeros, ICoreComplexMatrix& outPoles,
double& outGain, std::string* whyNot = nullptr);
// MATLAB's [bz, az] = impinvar(b, a, fs[, tol]): the discrete filter whose
// impulse response is the analog one SAMPLED at fs -- a different
// transformation from the bilinear one and not a variant of it, since it
// aliases rather than warping.
//
// The route is `impinvar.m`'s, and its shape is the reason the answer is
// what it is: the analog filter is expanded into partial fractions, each
// pole p contributes exp(p/fs) to the discrete denominator, the first
// `order` samples of the analog impulse response are formed directly, and
// the numerator is read back out by filtering that response through the
// new denominator. A repeated pole contributes t^(m-1)/(m-1)! alongside
// its exponential, which is where `tol` comes in: it is the RELATIVE
// distance within which two poles are one repeated pole (MATLAB's 1e-3).
static bool impulseInvariant(const ICoreMatrix& numerator, const ICoreMatrix& denominator,
const double& sampleRate, const double& tolerance,
ICoreMatrix& outNumerator, ICoreMatrix& outDenominator,
std::string* whyNot = nullptr);
};
};ICoreFirDesign.h#
src/ICoreBlocks/ICoreMath/SignalProcessing/ICoreFirDesign.h
ICoreFirDesign#
ICoreFirDesign.h:38 · class · nested OrderEstimate · 0 declaration(s)
The FIR half of filter design (toolbox row T2.4): firls, fir1, fir2, kaiserord and firpmord.
class ICoreFirDesign {
public:
// Which linear-phase family `firls` is asked for. `Standard` is the
// symmetric Type I/II filter every ordinary design wants.
enum class LeastSquaresType { Standard, Hilbert, Differentiator };
// MATLAB's b = firls(n, f, a[, w][, ftype]). `frequencies` runs 0..1 in
// PAIRS (Nyquist is 1), `amplitudes` gives the template's value at each
// edge, and `weights` holds one weight per band (empty for all ones).
//
// The answer is `order + 1` taps, or one more when firchk's rule above
// applies -- which is reported through `orderRaised` rather than warned
// about, since this console cannot warn from inside an expression.
static bool leastSquares(const size_t& order, const ICoreMatrix& frequencies,
const ICoreMatrix& amplitudes, const ICoreMatrix* weights,
const LeastSquaresType& type, ICoreMatrix& taps,
bool* orderRaised = nullptr, std::string* whyNot = nullptr);
// Which band arrangement `fir1` is asked for. `Low` and `Bandpass` are the
// defaults for a scalar and a two-entry Wn; `High` and `Stop` are the
// words, and a longer Wn is a multiband whose first band follows the word.
enum class Band { Low, High, Bandpass, Stop, DcZero, DcOne };
// MATLAB's b = fir1(n, Wn[, ftype][, window][, 'scale' | 'noscale']).
// `Wn` is normalised to NYQUIST (1 is half the sample rate), `window` is
// the n+1 window to multiply by (empty asks for MATLAB's Hamming default),
// and `scale` is MATLAB's default: the filter is divided by its own gain
// at the centre of the first band, so a lowpass has unit gain at DC.
static bool windowed(const size_t& order, const ICoreMatrix& cutoffs, const Band& band,
const ICoreMatrix* window, const bool& scale, ICoreMatrix& taps,
std::string* whyNot = nullptr);
// MATLAB's b = fir2(n, f, m[, npt][, lap][, window]): the frequency
// SAMPLING design -- interpolate the template onto a dense grid, give it
// the linear phase of an n-tap filter, transform back and window.
//
// `points` of 0 asks for MATLAB's default (512 below 1024 taps, the next
// power of two above that); `ramp` of 0 asks for its `fix(npt / 25)`,
// which is how wide a jump in the template is smeared before transforming.
static bool frequencySampled(const size_t& order, const ICoreMatrix& frequencies,
const ICoreMatrix& magnitudes, const size_t& points,
const size_t& ramp, const ICoreMatrix* window,
ICoreMatrix& taps, std::string* whyNot = nullptr);
// What an order estimate answers: the order, the band edges it wants a
// design to use, and (for kaiserord) the window parameter beta.
struct OrderEstimate {
size_t order = 0;
ICoreMatrix cutoffs; // normalised to Nyquist, ready for fir1
double beta = 0.0; // kaiserord only
Band band = Band::Low; // which ftype the edges imply
ICoreMatrix frequencies; // firpmord only: the firpm band vector
ICoreMatrix amplitudes; // firpmord only
ICoreMatrix weights; // firpmord only
};
// MATLAB's [n, Wn, beta, ftype] = kaiserord(f, a, dev[, fs]): Kaiser's own
// closed-form estimate of how long a windowed design must be to hold a
// given ripple across a given transition. `sampleRate` of 0 means the
// edges are already normalised to a rate of 2 (MATLAB's default).
static bool kaiserOrder(const ICoreMatrix& edges, const ICoreMatrix& amplitudes,
const ICoreMatrix& deviations, const double& sampleRate,
OrderEstimate& out, std::string* whyNot = nullptr);
// MATLAB's b = firpm(n, f, a[, w]): the EQUIRIPPLE design -- the filter
// whose largest weighted error is as small as it can be, which is a
// different filter from firls's (least squares) and usually a shorter one
// for the same specification.
//
// It is the only iterative design in this file. The Remez exchange guesses
// `nfcns + 1` extremal frequencies, fits the unique response that
// alternates through them, then moves each one to the nearby grid point
// where the error is actually largest, and repeats. What makes it
// reproducible rather than approximate is that the exchange is DISCRETE:
// the extrema only ever move to points of a fixed frequency grid, so two
// implementations that build the same grid and apply the same moves agree
// to rounding rather than to a tolerance.
//
// Three things about it are transcribed rather than derived, because none
// of them follows from the mathematics:
//
// THE GRID IS PART OF THE ANSWER. `firpmgrid` lays out `1/(16*nfcns)`
// spacing per band, forces the last point of each band onto the band
// edge, and drops or moves the final point when the last band reaches
// Nyquist. A denser grid is a different filter, not a better-converged
// one, and MATLAB quadruples the density only when the grid came out
// shorter than the filter.
//
// THE BARYCENTRIC WEIGHTS SKIP EVERY `jet`-th POINT (`remezdd`), where
// `jet = (nfcns - 1)/15 + 1`. That is a numerical-conditioning trick of
// the original FORTRAN, not a definition, and it changes the last bits.
//
// THE COEFFICIENTS COME OUT THROUGH A CHEBYSHEV CHANGE OF VARIABLE that
// is SKIPPED when the design spans the whole band or the filter is very
// short (`kkk`), and taken when it does not. Both branches are here
// because MATLAB takes both.
//
// `deviation` is the converged ripple -- MATLAB's second output -- and is
// the number that says whether the specification was met at this order.
// The Hilbert and differentiator types are refused, as `leastSquares`
// refuses them, and for the same reason.
static bool equiripple(const size_t& order, const ICoreMatrix& frequencies,
const ICoreMatrix& amplitudes, const ICoreMatrix* weights,
const LeastSquaresType& type, ICoreMatrix& taps,
double* deviation = nullptr, std::string* whyNot = nullptr);
// MATLAB's [B, G] = sgolay(order, frameLength[, weights]): the
// Savitzky-Golay projection and its differentiation matrix.
//
// B is the projector onto the frame's polynomials -- row (frameLength+1)/2
// smooths the middle of a frame and the rows either side of it smooth each
// edge position -- and it is `ICoreSignalFilters`'s matrix rather than a
// second copy, because `sgolayfilt` is nothing without it and two copies
// would eventually disagree.
//
// G is S * inv(S' * S) where S is the frame's Vandermonde matrix. MATLAB
// spells it Q/R' off an economy QR, which is the SAME matrix by a
// different route -- so it does not depend on which QR produced it, and
// this computes it from the normal equations instead. Column k+1 of G is
// the least-squares estimate of the k-th derivative (times k!), which is
// what makes the second output worth having at all.
static bool savitzkyGolay(const size_t& order, const size_t& frameLength,
const ICoreMatrix* weights, ICoreMatrix& projection,
ICoreMatrix* differentiation, std::string* whyNot = nullptr);
// MATLAB's [n, fo, ao, w] = firpmord(f, a, dev[, fs]): the same question
// for an EQUIRIPPLE design, by Herrmann's formula rather than Kaiser's --
// a different estimate of a different filter, and the two do not agree.
static bool remezOrder(const ICoreMatrix& edges, const ICoreMatrix& amplitudes,
const ICoreMatrix& deviations, const double& sampleRate,
OrderEstimate& out, std::string* whyNot = nullptr);
};
};ICoreFrequencyResponse.h#
src/ICoreBlocks/ICoreMath/SignalProcessing/ICoreFrequencyResponse.h
Both are held BY VALUE in Response below, so a forward declaration will not do: a caller of this header needs their size.
ICoreFrequencyResponse#
ICoreFrequencyResponse.h:39 · class · nested Response · 0 declaration(s)
The frequency response of a filter, digital and analog (toolbox row T2.6): what a design is LOOKED AT with, once the prototypes and the transformations have produced one.
class ICoreFrequencyResponse {
public:
// The response and the grid it was taken on, which travel together because
// a response without its frequencies is a picture nobody can place.
struct Response {
ICoreComplexMatrix values;
ICoreMatrix frequencies;
};
// MATLAB's freqz on a coefficient pair. Exactly one of `frequencies` and
// `points` is used: a non-null `frequencies` is the caller's own grid
// (returned unchanged), and otherwise `points` picks one of the four
// spellings above -- 0 asks for MATLAB's default of 512.
//
// `sampleRate` of 0 asks for the normalized grid in radians per sample;
// any positive value answers in Hz. `whole` swaps 0, pi) for the whole
// circle.
static bool digital(const ICoreMatrix& numerator, const ICoreMatrix& denominator,
const ICoreMatrix* frequencies, const size_t& points,
const bool& whole, const double& sampleRate,
Response& out, std::string* whyNot = nullptr);
// MATLAB's freqz(sos, ...): the same response for a filter given as an
// L x 6 matrix of second-order sections, whose response is the PRODUCT of
// the sections' -- so a design that was split into sections for numerical
// reasons is looked at without being reassembled into one long pair.
static bool sections(const ICoreMatrix& sos, const ICoreMatrix* frequencies,
const size_t& points, const bool& whole, const double& sampleRate,
Response& out, std::string* whyNot = nullptr);
// MATLAB's freqs(b, a, w): the ANALOG response, b(jw)/a(jw). The frequency
// vector is not optional here and not a count -- an analog filter has no
// Nyquist rate to spread a default grid over.
static bool analog(const ICoreMatrix& numerator, const ICoreMatrix& denominator,
const ICoreMatrix& frequencies, ICoreComplexMatrix& out,
std::string* whyNot = nullptr);
};
};ICoreIirDesign.h#
[src/ICoreBlocks/ICoreMath/SignalProcessing/ICoreIirDesign.h
ICoreIirDesign#
ICoreIirDesign.h:53 · class · nested Design · 0 declaration(s)
The five classical IIR designs and the four order estimates that size them (toolbox rows T2.1, T2.2 and T2.3): butter, cheby1, cheby2, ellip, besself, and buttord, cheb1ord, `cheb2ord...
class ICoreIirDesign {
public:
// Which approximation, which is to say which prototype the pipeline
// starts from. `Bessel` is analog-only, as MATLAB's `besself` is.
enum class Approximation { Butterworth, ChebyshevI, ChebyshevII, Elliptic, Bessel };
// MATLAB's own btype numbering, which the band edges also decide: one
// edge is a lowpass unless `high` says otherwise, two edges are a
// bandpass unless `stop` does.
enum class Band { Lowpass = 1, Bandpass = 2, Highpass = 3, Bandstop = 4 };
// A designed filter in the form every one of MATLAB's five answers it in.
// The zeros are as many as the poles for a discrete design (the bilinear
// map sends infinity to z = -1) and fewer for an analog lowpass.
struct Design {
std::vector<std::complex<double>> zeros;
std::vector<std::complex<double>> poles;
double gain = 1.0;
};
// The five designs. `edges` is MATLAB's Wn / Wp / Ws / Wo: one entry for a
// lowpass or highpass, two for a bandpass or bandstop, normalised to
// Nyquist (strictly inside 0..1) unless `analog`, where it is rad/s.
//
// `passbandRipple` is read by ChebyshevI and Elliptic, `stopbandAttenuation`
// by ChebyshevII and Elliptic; the other approximations ignore both, which
// is what lets the five share one signature.
static bool design(const Approximation& approximation, const size_t& order,
const double& passbandRipple, const double& stopbandAttenuation,
const std::vector<double>& edges, const Band& band, const bool& analog,
Design& out, std::string* whyNot = nullptr);
// MATLAB's [b, a] out of a design: `den = real(poly(p))` and
// `num = [zeros(1, np - nz) k*real(poly(z))]`, both in descending powers
// and both of length np + 1 -- the numerator is LEADING-ZERO PADDED rather
// than shortened, so `b` and `a` are always the same length and a filter
// of order n reads as one.
static void coefficients(const Design& design, std::vector<double>& numerator,
std::vector<double>& denominator);
// The four order estimates: the least order whose filter is inside `rp` dB
// of ripple across the passband and at least `rs` dB down across the
// stopband, and the natural frequency (or pair) to hand the design.
//
// Each approximation reads its own natural frequency off a different
// place, which is MATLAB's own arrangement and not a tidiness this could
// fix: `buttord` answers the -3 dB frequency it computed, `cheb1ord` the
// PASSBAND edges it was given, `cheb2ord` the STOPBAND edges it was given,
// and `ellipord` the passband edges again. `Bessel` has no estimate --
// MATLAB has no `besselord` -- and is refused by name.
static bool orderEstimate(const Approximation& approximation,
const std::vector<double>& passbandEdges,
const std::vector<double>& stopbandEdges,
const double& passbandRipple, const double& stopbandAttenuation,
const bool& analog, size_t& order,
std::vector<double>& naturalEdges,
std::string* whyNot = nullptr);
};
};ICoreLinearPrediction.h#
src/ICoreBlocks/ICoreMath/SignalProcessing/ICoreLinearPrediction.h
ICoreLinearPrediction#
ICoreLinearPrediction.h:41 · class · 0 declaration(s)
Linear prediction: the Signal Processing Toolbox family that turns a signal (or its autocorrelation) into an all-pole model, and converts between the three equivalent spellings of one (toolbox row ...
class ICoreLinearPrediction {
public:
// MATLAB's [a, e, k] = levinson(r, order): the Levinson-Durbin recursion
// over the autocorrelation sequence `r` (r(0) first, whichever way the
// vector is turned).
//
// `a` comes back as a ROW of order+1 coefficients beginning with 1, `error`
// is the FINAL prediction error E_p, and `reflection` is a COLUMN of the
// order reflection coefficients -- MATLAB's three shapes exactly.
//
// `order` may not exceed numel(r) - 1; MATLAB clamps it to that, and so
// does this, because a recursion has nothing left to read past the end of
// its own data.
static bool levinson(const ICoreMatrix& r, const size_t& order,
ICoreMatrix& a, double& error, ICoreMatrix& reflection,
std::string* whyNot = nullptr);
// MATLAB's [k, r0] = poly2rc(a, efinal): the step-DOWN recursion, from the
// prediction-error polynomial to the reflection coefficients. `efinal` is
// the final prediction error and `zeroLagError` comes back as the
// zero-lag error r0 = efinal / prod(1 - k_i^2).
//
// A reflection coefficient of magnitude exactly 1 stops the recursion --
// its step divides by 1 - k^2 -- and is refused with that reason rather
// than answered as an infinity.
static bool polyToReflection(const ICoreMatrix& a, const double& finalError,
ICoreMatrix& reflection, double& zeroLagError,
std::string* whyNot = nullptr);
// MATLAB's [a, efinal] = rc2poly(k, r0): the step-UP recursion, the exact
// inverse of the one above. `r0` is the ZERO-LAG error and `finalError`
// comes back as r0 * prod(1 - k_i^2).
static bool reflectionToPoly(const ICoreMatrix& reflection, const double& zeroLagError,
ICoreMatrix& a, double& finalError,
std::string* whyNot = nullptr);
// MATLAB's r = poly2ac(a, efinal): the autocorrelation sequence the
// polynomial came from, as a COLUMN of numel(a) entries beginning with r0.
// Its inverse is `levinson` at full order, which is MATLAB's `ac2poly`.
static bool polyToAutocorrelation(const ICoreMatrix& a, const double& finalError,
ICoreMatrix& r, std::string* whyNot = nullptr);
// MATLAB's lpc(x, p) and aryule(x, p), which are one estimator: the
// Levinson recursion over the BIASED autocorrelation of `x` (divided by
// numel(x), so the model is windowed and always stable). `a` is a row,
// `error` the final prediction error -- lpc's `g` and aryule's `e` are the
// same number -- and `reflection` a column, which only aryule reports.
static bool yuleWalker(const ICoreMatrix& x, const size_t& order,
ICoreMatrix& a, double& error, ICoreMatrix& reflection,
std::string* whyNot = nullptr);
// MATLAB's [a, e, k] = arburg(x, p): Burg's recursion on the forward and
// backward prediction errors. Same three shapes as levinson's.
static bool burg(const ICoreMatrix& x, const size_t& order,
ICoreMatrix& a, double& error, ICoreMatrix& reflection,
std::string* whyNot = nullptr);
// MATLAB's arcov(x, p) and, with `modified` set, armcov(x, p): the
// least-squares solve over the covariance data matrix -- the forward
// errors alone, or the forward and backward ones together.
//
// `error` is the white-noise variance estimate, which is a MEAN square:
// the data matrix carries MATLAB's own 1/sqrt(N - p) scaling (and a
// further 1/sqrt(2) for the modified method), and dropping it would answer
// a model with the right coefficients and the wrong variance.
static bool covarianceMethod(const ICoreMatrix& x, const size_t& order,
const bool& modified, ICoreMatrix& a, double& error,
std::string* whyNot = nullptr);
};
};ICoreParametricSpectrum.h#
src/ICoreBlocks/ICoreMath/SignalProcessing/ICoreParametricSpectrum.h
ICoreParametricSpectrum#
ICoreParametricSpectrum.h:41 · class · nested Spectrum · 0 declaration(s)
The parametric (autoregressive) spectral estimates: pyulear, pburg, pcov and pmcov (toolbox row T2.18).
class ICoreParametricSpectrum {
public:
// Which all-pole fit the spectrum is taken of. The four are MATLAB's four
// names, and the difference between them is entirely in the fit.
enum class Estimator {
YuleWalker, // pyulear -- aryule
Burg, // pburg -- arburg
Covariance, // pcov -- arcov
ModifiedCovariance, // pmcov -- armcov
};
// Which part of the circle is reported. `OneSided` is MATLAB's default for
// a real signal and folds as described above; `TwoSided` reports the whole
// circle unfolded.
enum class Range { OneSided, TwoSided };
// The density and the grid it was taken on, which travel together for the
// same reason a frequency response and its grid do.
struct Spectrum {
ICoreMatrix values;
ICoreMatrix frequencies;
};
// MATLAB's [pxx, w] = pyulear/pburg/pcov/pmcov(x, order[, nfft][, fs]).
// `points` of 0 asks for MATLAB's default of 256; `sampleRate` of 0 asks
// for the normalised grid in radians per sample, and any positive value
// answers in Hz with the density scaled by the rate instead of by 2*pi.
//
// An order at or above the signal's length is refused, as MATLAB refuses
// it: there is no model left to fit.
static bool estimate(const ICoreMatrix& signal, const size_t& order, const size_t& points,
const Estimator& estimator, const Range& range,
const double& sampleRate, Spectrum& out,
std::string* whyNot = nullptr);
};
};ICorePeakFinding.h#
src/ICoreBlocks/ICoreMath/SignalProcessing/ICorePeakFinding.h
ICorePeakFinding#
ICorePeakFinding.h:43 · class · nested Options, Peaks · 0 declaration(s)
MATLAB's findpeaks (toolbox row T2.23): which samples of a signal are local maxima, how far each one stands above the terrain around it, and how wide it is at half that height.
class ICorePeakFinding {
public:
// The reference level a peak's width is measured at.
enum class WidthReference {
HalfProminence, // MATLAB's default 'halfprom': (peak + base) / 2
HalfHeight, // 'halfheight': peak / 2, bounded by the saddles
};
// What order the surviving peaks are reported in. `None` is MATLAB's
// default and reports them in the order they occur.
enum class Sort { None, Ascend, Descend };
// Every thinning rule, with MATLAB's defaults -- so a default-constructed
// Options is `findpeaks(y)`.
struct Options {
double minHeight = -std::numeric_limits<double>::infinity();
double threshold = 0.0;
double minProminence = 0.0;
double minWidth = 0.0;
double maxWidth = std::numeric_limits<double>::infinity();
double minDistance = 0.0; // measured in the LOCATION units, not samples
size_t maxPeaks = 0; // 0 asks for every peak that survives
Sort sort = Sort::None;
WidthReference widthReference = WidthReference::HalfProminence;
};
// The four outputs of `[pks, locs, w, p] = findpeaks(...)`, which travel
// together because three of them are only meaningful beside the first.
// Each is shaped the way `y` was -- a row signal answers rows.
struct Peaks {
ICoreMatrix values;
ICoreMatrix locations;
ICoreMatrix widths;
ICoreMatrix prominences;
};
// MATLAB's findpeaks on a real vector. `locations` is the x axis the peaks
// are reported and the minimum distance is measured on: pass a null
// pointer for MATLAB's default of 1..n (which is what `findpeaks(y)` uses)
// or a strictly increasing vector of the same length (`findpeaks(y, x)`;
// a sample rate is (0:n-1)/fs, formed by the caller).
//
// A signal shorter than three samples is refused, as MATLAB refuses it:
// the first and last sample can never be strictly greater than a neighbour
// on both sides, so there is nothing for the definition to answer.
static bool find(const ICoreMatrix& signal, const ICoreMatrix* locations,
const Options& options, Peaks& out, std::string* whyNot = nullptr);
};
};ICoreResampling.h#
src/ICoreBlocks/ICoreMath/SignalProcessing/ICoreResampling.h
ICoreResampling#
ICoreResampling.h:41 · class · 0 declaration(s)
The rate changers (toolbox row T2.21): upsample, downsample, upfirdn, interp, decimate and resample.
class ICoreResampling {
public:
// MATLAB's y = upsample(x, n[, phase]) and y = downsample(x, n[, phase]):
// exact, and each other's inverse when the phase agrees. The answer keeps
// the input's orientation, as MATLAB's does.
static bool upsample(const ICoreMatrix& x, const size_t& factor, const size_t& phase,
ICoreMatrix& out, std::string* whyNot = nullptr);
static bool downsample(const ICoreMatrix& x, const size_t& factor, const size_t& phase,
ICoreMatrix& out, std::string* whyNot = nullptr);
// MATLAB's y = upfirdn(x, h, p, q): x upsampled by `p`, filtered by `h`,
// then decimated by `q`, in one pass.
//
// **The output LENGTH is the part that is not obvious** and it is the
// trap the board's note names: `ceil(((Lx - 1)*p + Lh) / q)`, the length of
// the full convolution of the upsampled signal, thinned by q -- so a
// filter longer than the signal still lengthens the answer, and the last
// output sample can be one the filter has not finished.
static bool upfirdn(const ICoreMatrix& x, const ICoreMatrix& h, const size_t& p,
const size_t& q, ICoreMatrix& out, std::string* whyNot = nullptr);
// MATLAB's y = interp(x, r[, n[, cutoff]]): the signal interpolated to r
// times as many samples, so that y(1:r:end) IS x -- which the FIR design
// below guarantees and a general lowpass does not.
//
// `n` is the number of ORIGINAL samples on each side that each new sample
// is computed from (MATLAB's default 4, so a 2*4*r+1 tap filter) and
// `cutoff` the normalised bandwidth (0.5). The signal must be longer than
// 2*n samples, which is MATLAB's own condition and is refused with it.
static bool interpolate(const ICoreMatrix& x, const size_t& factor, const size_t& neighbours,
const double& cutoff, ICoreMatrix& out,
std::string* whyNot = nullptr);
// MATLAB's y = decimate(x, r[, n][, "fir"]): every r-th sample of the
// low-passed signal, `ceil(numel(x)/r)` of them.
//
// The two arms are different filters AND different edge treatments, which
// is why they answer different numbers: the Chebyshev I arm is zero-phase
// (`filtfilt`, so no delay to correct) and the FIR arm is a single forward
// pass whose group delay is corrected by starting the read at
// `round(gd(1) + 1.25)`. ⚠ MATLAB reduces the IIR order until the filter
// it designs really is `rip` dB down at 0.8/r -- a high order at a large r
// is numerically dead, and the loop is reproduced rather than the order
// being trusted.
static bool decimate(const ICoreMatrix& x, const size_t& factor, const size_t& order,
const bool& finiteImpulse, ICoreMatrix& out,
std::string* whyNot = nullptr);
// MATLAB's y = resample(x, p, q[, n[, beta]]): `ceil(numel(x)*p/q)`
// samples of x at p/q times its rate, through a Kaiser-windowed
// least-squares FIR.
//
// `p` and `q` are reduced to lowest terms first, as MATLAB reduces them
// (`rat(p/q, 1e-12)`), so resample(x, 4, 6) IS resample(x, 2, 3) on both
// sides. `n` of 0 asks for the degenerate all-ones filter MATLAB uses
// there; `outFilter`, when not null, receives the filter actually used --
// MATLAB's second output.
static bool resample(const ICoreMatrix& x, const size_t& up, const size_t& down,
const size_t& neighbours, const double& beta, ICoreMatrix& out,
ICoreMatrix* outFilter, std::string* whyNot = nullptr);
};
};ICoreSignalFilters.h#
src/ICoreBlocks/ICoreMath/SignalProcessing/ICoreSignalFilters.h
ICoreSignalFilters#
ICoreSignalFilters.h:38 · class · 0 declaration(s)
The five Signal Processing filtering helpers of toolbox row T2.11 -- sosfilt, fftfilt, medfilt1, sgolayfilt and hampel.
class ICoreSignalFilters {
public:
// MATLAB's `sosfilt(sos, x)`: a cascade of biquads, one per ROW of the
// L x 6 `sos` matrix ([b0 b1 b2 a0 a1 a2]). Each row is divided through by
// its own a0 and each stage runs over the previous stage's output.
//
// There is no gain vector here -- that is `filtfilt(sos, g, x)`'s and
// `tf2sos`'s (T2.8) argument, not this one. `x` is filtered down its
// columns, or along itself when it is a vector, and the answer has x's
// shape.
static ICoreMatrix sections(const ICoreMatrix& sos, const ICoreMatrix& x,
std::string* whyNot = nullptr);
// MATLAB's `fftfilt(b, x)` and `fftfilt(b, x, nfft)`: the FIR filter `b`
// applied by overlap-add, answering the FIRST numel(x) samples of the
// convolution -- the same numbers `filter(b, 1, x)` answers, reached by a
// different route and therefore not to the last bit.
//
// `nfft <= 0` means "choose it", and the choice is MATLAB's own: the
// published flop table in fftfilt.m, minimising ceil(nx/L)*flops(nfft)
// over the power-of-two lengths that leave a usable block. A given `nfft`
// is raised to at least numel(b) and then to the next power of two, as
// MATLAB raises it.
static ICoreMatrix blockConvolution(const ICoreMatrix& b, const ICoreMatrix& x,
const long long& nfft, std::string* whyNot = nullptr);
// How `medfilt1` treats the samples a window reaches past the ends.
// MATLAB's default is ZeroPad -- the window stays full and the missing
// samples are zeros, which pulls the answer towards zero at both ends.
// Truncate is `movmedian`'s (and this console's) shrinking window.
enum class Endpoints { ZeroPad, Truncate };
// MATLAB's `medfilt1(x, n)`: the running median of an n-wide window.
// An even `n` reaches one further BACK than forward (floor(n/2) back,
// ceil(n/2)-1 forward), which is `movmedian`'s rule and therefore the
// console's. `n == 0` answers all-NaN, as medfilt1.m does.
//
// `omitNaN` is MATLAB's 'omitnan'; the default there and here is
// 'includenan', under which one NaN anywhere in a window makes that
// window's answer NaN.
static ICoreMatrix runningMedian(const ICoreMatrix& x, const size_t& n,
const Endpoints& ends, const bool& omitNaN,
std::string* whyNot = nullptr);
// `sgolay`'s first output B (T2.4 owns the NAME; this is the matrix, kept
// here because sgolayfilt is nothing without it and two copies of it would
// eventually disagree). B is frameLength x frameLength: row m+1 is the
// smoothing filter for the middle of the frame, and the rows either side
// of it are the filters for each edge position.
//
// `weights` is NULL for unity. It is a pointer rather than an empty matrix
// because a default-constructed ICoreMatrix is 1 x 1 holding 0 and not
// empty -- numel() cannot tell "no weights" from "one weight", and reading
// it as the latter is what made every unweighted sgolayfilt call refuse
// itself. It is the least-squares weighting, and MATLAB folds it in as
// sqrt(w) on both sides of the projector.
static bool savitzkyGolayProjection(const size_t& order, const size_t& frameLength,
const ICoreMatrix* weights, ICoreMatrix& B,
std::string* whyNot = nullptr);
// MATLAB's `sgolayfilt(x, order, frameLength[, weights])`. `frameLength`
// must be ODD and greater than `order`, and x must be at least
// frameLength long -- MATLAB's own three refusals.
//
// The middle of the answer is one FIR pass with row (frameLength+1)/2 of
// B; the first and last (frameLength-1)/2 samples are not that filter
// started early, they are the polynomial fitted to the first (last) whole
// frame, evaluated at each edge position.
static ICoreMatrix savitzkyGolayFilter(const ICoreMatrix& x, const size_t& order,
const size_t& frameLength,
const ICoreMatrix* weights,
std::string* whyNot = nullptr);
// MATLAB's `[y, i, xmedian, xsigma] = hampel(x, k, nsigma)`: the outlier
// filter that replaces a sample by its local median when it is more than
// `nsigma` local standard deviations away from it.
//
// The window is 2k+1 wide with SHRINKING endpoints, `xsigma` is the moving
// MEDIAN absolute deviation scaled by 1/(sqrt(2)*erfcinv(3/2)) = 1.4826...
// (so it estimates sigma for Gaussian data), and `outliers` is MATLAB's
// logical mask as 1s and 0s. The comparison is written NEGATED -- a sample
// is an outlier when NOT(|x - med| <= nsigma*sigma) -- so a NaN sample is
// an outlier, which a straight `>` would answer the other way.
static bool outlierFilter(const ICoreMatrix& x, const size_t& k, const double& nsigma,
ICoreMatrix& filtered, ICoreMatrix& outliers,
ICoreMatrix& median, ICoreMatrix& sigma,
std::string* whyNot = nullptr);
};
};ICoreSignalQuality.h#
src/ICoreBlocks/ICoreMath/SignalProcessing/ICoreSignalQuality.h
ICoreSignalQuality#
ICoreSignalQuality.h:45 · class · nested Measurement · 0 declaration(s)
The four single-tone quality measurements (the toolbox board's T2.25): snr, thd, sinad and sfdr -- how much of a signal's power is the tone you meant, and how much is everything else.
class ICoreSignalQuality {
public:
// What a measurement answers. Not every field is filled by every name --
// each method's comment says which -- and the ratio is always in dB.
struct Measurement {
double ratio = 0.0; // r, in dB
double power = 0.0; // the noise/distortion/spur power, dB
double frequency = 0.0; // sfdr's spur frequency
std::vector<double> harmonicPower; // thd's per-harmonic power, dB
std::vector<double> harmonicFrequency; // and where each one was found
};
// snr(x, fs[, n]) -- the fundamental's power over everything that is not
// the fundamental, its first `n` harmonics (default 6) or DC. Fills
// `ratio` and `power` (the total noise, dB).
//
// `aliased` is MATLAB's 'aliased' word: a harmonic above Nyquist is folded
// back onto the band rather than ignored.
static bool signalToNoise(const ICoreMatrix& signal, const double& sampleRate,
const size_t& harmonics, const bool& aliased,
Measurement& out, std::string* whyNot);
// snr(x, y) -- the two-signal form, which has no spectrum in it at all:
// 10*log10(sum(x.^2) / sum(y.^2)). MATLAB tells this form from the one
// above by y being a VECTOR rather than a scalar sample rate, and so does
// the console.
static bool signalToNoise(const ICoreMatrix& signal, const ICoreMatrix& noise,
Measurement& out, std::string* whyNot);
// thd(x, fs[, n]) -- the summed power of harmonics 2..n over the
// fundamental's. Fills `ratio`, and `harmonicPower`/`harmonicFrequency`
// with the fundamental FIRST (MATLAB's order), each power in dB.
static bool totalHarmonicDistortion(const ICoreMatrix& signal, const double& sampleRate,
const size_t& harmonics, const bool& aliased,
Measurement& out, std::string* whyNot);
// sinad(x, fs) -- the fundamental over everything else, harmonics
// INCLUDED, which is the one line that separates it from snr: it removes
// DC and the fundamental and calls the whole remainder noise. Fills
// `ratio` and `power`.
static bool signalToNoiseAndDistortion(const ICoreMatrix& signal, const double& sampleRate,
Measurement& out, std::string* whyNot);
// sfdr(x, fs[, msd]) -- the fundamental over the largest SPUR, whatever it
// is and wherever it is; `msd` is a minimum spurious distance, a band
// around the fundamental to exclude. Fills `ratio`, `power` (the spur's,
// dB) and `frequency`.
static bool spuriousFreeDynamicRange(const ICoreMatrix& signal, const double& sampleRate,
const double& minimumSpurDistance, Measurement& out,
std::string* whyNot);
};
};ICoreSpectralEstimation.h#
src/ICoreBlocks/ICoreMath/SignalProcessing/ICoreSpectralEstimation.h
ICoreSpectralEstimation#
ICoreSpectralEstimation.h:44 · class · nested Options, Spectrum, ShortTime · 0 declaration(s)
The non-parametric spectral estimates (toolbox rows T2.14–T2.17): periodogram, pwelch, cpsd, mscohere, tfestimate and spectrogram.
class ICoreSpectralEstimation {
public:
// Whether the answer is a DENSITY (power per unit frequency, MATLAB's
// 'psd') or a POWER (per bin, 'power'). The difference is the normaliser.
enum class Scaling { Density, Power };
// Which part of the circle is reported. `OneSided` folds and is MATLAB's
// default for a real signal; `Centered` puts DC in the middle.
enum class Range { OneSided, TwoSided, Centered };
// What the six names ask this machine for. `Auto` is one signal;
// `Cross`, `Coherence` and `Transfer` take two.
enum class Estimate { Auto, Cross, Coherence, Transfer };
struct Options {
size_t points = 0; // nfft; 0 asks for the estimate's default
double sampleRate = 0.0; // 0 answers in radians per sample (Fs = 2*pi)
Range range = Range::OneSided;
Scaling scaling = Scaling::Density;
};
// The estimate and the grid it was taken on, which travel together.
//
// Two of the four estimates are COMPLEX -- a cross spectrum's phase is the
// lag between the two signals and a transfer estimate's is the system's,
// which is most of what either is asked for -- so the answer carries both
// shapes and says which one it filled. A coherence and an auto spectrum
// are real.
struct Spectrum {
ICoreMatrix values; // filled when isComplex is false
ICoreComplexMatrix complexValues; // filled when isComplex is true
bool isComplex = false;
ICoreMatrix frequencies;
};
// MATLAB's [pxx, f] = periodogram(x[, window[, nfft[, fs]]]): one segment,
// one transform. `window` may be null for MATLAB's default rectangular
// window of the signal's own length; `points` of 0 is
// `max(256, 2^nextpow2(n))`.
static bool periodogram(const ICoreMatrix& signal, const ICoreMatrix* window,
const Options& options, Spectrum& out,
std::string* whyNot = nullptr);
// MATLAB's pwelch / cpsd / mscohere / tfestimate: the same estimate over
// overlapping segments, averaged.
//
// `window` is the caller's window vector, or null; `windowLength` is a
// scalar window argument (0 for none), which asks for a Hamming window of
// that length. `overlap` is used only when `overlapGiven` is set.
//
// `second` is the other signal for the three two-signal estimates and is
// null for `Auto`. The ORDER matters and is MATLAB's: a transfer estimate
// is the cross spectrum of (y, x) over the auto spectrum of x, so a
// reversed pair answers its reciprocal rather than an error.
static bool welch(const ICoreMatrix& signal, const ICoreMatrix* second,
const ICoreMatrix* window, const size_t& windowLength,
const bool& overlapGiven, const size_t& overlap,
const Estimate& estimate, const Options& options,
Spectrum& out, std::string* whyNot = nullptr);
// The short-time transform behind `spectrogram`: one column per segment,
// the segments' centre times, and the frequency grid.
//
// `values` is the complex STFT — MATLAB's first output — and `power` is
// its PSD (the fourth output), which is not `abs(values)^2`: it carries
// the same normalisation and fold the estimates above do.
struct ShortTime {
ICoreComplexMatrix values;
ICoreMatrix frequencies;
ICoreMatrix times;
ICoreMatrix power;
};
static bool shortTime(const ICoreMatrix& signal, const ICoreMatrix* window,
const size_t& windowLength, const bool& overlapGiven,
const size_t& overlap, const Options& options,
ShortTime& out, std::string* whyNot = nullptr);
};
};ICoreSpectralMeasures.h#
src/ICoreBlocks/ICoreMath/SignalProcessing/ICoreSpectralMeasures.h
ICoreSpectralMeasures#
ICoreSpectralMeasures.h:37 · class · nested Band · 0 declaration(s)
The measurements taken OF a signal or of its spectrum (toolbox row T2.24): rms, peak2peak, peak2rms, rssq, bandpower, meanfreq, medfreq, obw and powerbw.
class ICoreSpectralMeasures {
public:
// MATLAB's rms, peak2peak, peak2rms and rssq: exact arithmetic over the
// samples, no spectrum involved. `rms` is core MATLAB
// (toolbox/matlab/datafun/) rather than Signal Processing, and lands here
// because the core board is closed and has no row for it.
static bool rootMeanSquare(const ICoreMatrix& signal, double& out,
std::string* whyNot = nullptr);
static bool peakToPeak(const ICoreMatrix& signal, double& out,
std::string* whyNot = nullptr);
static bool peakToRms(const ICoreMatrix& signal, double& out,
std::string* whyNot = nullptr);
static bool rootSumSquare(const ICoreMatrix& signal, double& out,
std::string* whyNot = nullptr);
// MATLAB's bandpower(x) — the mean square, exactly, with no spectrum at
// all — and bandpower(x, fs, [flo fhi]), which integrates a Hamming
// periodogram over the band. `hasRange` picks between them.
static bool bandPower(const ICoreMatrix& signal, const double& sampleRate,
const bool& hasRange, const double& low, const double& high,
double& out, std::string* whyNot = nullptr);
// What the four spectrum measurements answer. Each takes the signal and a
// sample rate (0 for radians per sample).
struct Band {
double frequency = 0.0; // meanfreq / medfreq
double power = 0.0; // the power in the band each one measured
double bandwidth = 0.0; // obw / powerbw
double low = 0.0; // the band's edges, for obw and powerbw
double high = 0.0;
};
// MATLAB's meanfreq and medfreq: the spectrum's centre of mass, and the
// frequency that halves its power.
static bool meanFrequency(const ICoreMatrix& signal, const double& sampleRate,
Band& out, std::string* whyNot = nullptr);
static bool medianFrequency(const ICoreMatrix& signal, const double& sampleRate,
Band& out, std::string* whyNot = nullptr);
// MATLAB's obw(x[, fs[, [], P]]): the band holding P percent of the power
// (99 by default), and powerbw(x[, fs[, [], R]]): the band whose density
// stays within R dB of the peak.
//
// ⚠ **`drop` HAS NO DEFAULT HERE AND THE CALLER MUST PASS MATLAB'S**, which
// is `10*log10(2)` = 3.0102999566398120 dB and NOT 3. `powerbw.m:18` sets
// it as `R = 10*log10(1/2)` -- the EXACT half power -- while the phrase
// everyone repeats, including MATLAB's own documentation title, is "the
// 3-dB bandwidth". The difference between 10^(-3/10) = 0.501187 and 1/2 is
// 0.24% of the reference level and about 0.34% of the answer, which is
// small enough to look like rounding and far too large to be it: this
// console refused the name for an afternoon over exactly that gap, on the
// theory that the shipped function did not match its own published file.
// It matched; the default did not. An explicit R reproduces to the last
// digit at 3, 6 and 10 dB.
static bool occupiedBandwidth(const ICoreMatrix& signal, const double& sampleRate,
const double& percent, Band& out,
std::string* whyNot = nullptr);
static bool powerBandwidth(const ICoreMatrix& signal, const double& sampleRate,
const double& drop, Band& out, std::string* whyNot = nullptr);
};
};ICoreWaveforms.h#
src/ICoreBlocks/ICoreMath/SignalProcessing/ICoreWaveforms.h
ICoreWaveforms#
ICoreWaveforms.h:28 · class · 6 declaration(s)
The Signal Processing Toolbox waveform generators (T2.26): periodic waves, aperiodic pulses and the swept-frequency cosine.
class ICoreWaveforms {
public:
// sin(pi * x) and cos(pi * x), exact at the whole and half-whole numbers.
static double sinpi(const double& x);
static double cospi(const double& x);
// sawtooth(t) / sawtooth(t, width) -- the modified sawtooth: a rising ramp
// from -1 to 1 over the first `width` of each 2*pi period and a falling
// ramp back over the rest. width is in [0, 1] and defaults to 1; the
// triangle wave is width = 0.5. Refuses a width outside [0, 1].
static ICoreMatrix sawtooth(const ICoreMatrix& t, const double& width, std::string* whyNot);
// square(t) / square(t, duty) -- +1 for the first `duty` PERCENT of each
// 2*pi period and -1 for the rest; duty defaults to 50. MATLAB does not
// range-check duty (a duty of 200 is all +1), and neither does this.
static ICoreMatrix square(const ICoreMatrix& t, const double& duty);
// sinc(x) = sin(pi*x) / (pi*x), the NORMALISED sinc with sinc(0) = 1 and
// a zero at every non-zero whole number. Not sin(x)/x.
static ICoreMatrix sinc(const ICoreMatrix& x);
// rectpuls(t) / rectpuls(t, w) -- the unit-height rectangle of width `w`
// (default 1) centred on t = 0. MATLAB's own half-open convention: the
// pulse is 1 on [-w/2, w/2) and its left edge is included to within eps,
// which is why the edges are compared against eps rather than exactly.
static ICoreMatrix rectpuls(const ICoreMatrix& t, const double& width);
// tripuls(t) / tripuls(t, w) / tripuls(t, w, skew) -- the unit-height
// triangle of width `w` (default 1) centred on t = 0, its peak moved to
// skew * w/2. skew is in [-1, 1] and defaults to 0 (symmetric).
static ICoreMatrix tripuls(const ICoreMatrix& t, const double& width, const double& skew,
std::string* whyNot);
// gauspuls(t, fc, bw, bwr) -- the Gaussian-modulated sinusoidal RF pulse,
// `part` selecting which of MATLAB's three outputs is wanted:
// 0 yc, the in-phase pulse (the single-output answer)
// 1 ys, the quadrature pulse
// 2 ye, the envelope
// fc defaults to 1000 Hz, bw to 0.5 (fractional bandwidth), bwr to -6 dB.
static ICoreMatrix gauspuls(const ICoreMatrix& t, const double& centerHz,
const double& fractionalBandwidth, const double& bandwidthRefDb,
const int& part, std::string* whyNot);
// gauspuls('cutoff', fc, bw, bwr, tpr) -- the time at which the pulse
// envelope has fallen to `tpr` dB below its peak (tpr defaults to -60).
static bool gauspulsCutoff(const double& centerHz, const double& fractionalBandwidth,
const double& bandwidthRefDb, const double& trailingPulseDb,
double& cutoffOut, std::string* whyNot);
// chirp(t, f0, t1, f1, method, phi, quadraticType) -- the swept-frequency
// cosine that is f0 Hz at t = 0 and f1 Hz at t = t1. Defaults, MATLAB's:
// f0 = 0, t1 = 1, f1 = 100, phi = 0 degrees, method "linear".
// `method` is "linear", "quadratic" or "logarithmic"; `quadraticType` is
// "", "concave" or "convex" and is read only by the quadratic sweep,
// where it flips the sweep MATLAB would otherwise choose from the
// direction of f0 -> f1.
static ICoreMatrix chirp(const ICoreMatrix& t, const double& f0, const double& t1,
const double& f1, const std::string& method, const double& phaseDeg,
const std::string& quadraticType, std::string* whyNot);
// pulstran's SAMPLED-prototype half: the pulse train made by shifting a
// sampled prototype `p` (taken at `fs` samples per second) to each delay
// in `delays`, scaling by `amplitudes` when they are given, and summing.
// `method` is an interp1 method name. The CONTINUOUS half -- a function
// handle or a function name shifted per delay -- is not here: it must
// call back into the expression evaluator, which is above this layer, so
// the evaluator owns that arm and this class owns the arithmetic.
static ICoreMatrix pulstranSampled(const ICoreMatrix& t, const std::vector<double>& delays,
const std::vector<double>& amplitudes,
const ICoreMatrix& prototype, const double& sampleRate,
const std::string& method, std::string* whyNot);
};
};ICoreWindows.h#
src/ICoreBlocks/ICoreMath/SignalProcessing/ICoreWindows.h
ICoreWindows#
ICoreWindows.h:39 · class · 4 declaration(s)
The Signal Processing Toolbox spectral windows (T2.22): the sixteen names window dispatches to, each answering an N x 1 COLUMN.
class ICoreWindows {
public:
enum class Sampling { Symmetric, Periodic };
// The one entry point, so `window(@hann, n)` and `hann(n)` cannot drift:
// both come through here. `name` is the MATLAB spelling, lowercase.
//
// `parameter` is the second argument the four parameterised windows take
// and the other twelve refuse -- kaiser's beta (default 0.5), gausswin's
// alpha (2.5), tukeywin's ratio (0.5) and chebwin's sidelobe attenuation
// in dB (100). `hasParameter` says whether the caller gave one.
//
// `length` is a double because MATLAB takes one: a non-integer length is
// ROUNDED (MATLAB warns and rounds; this rounds), a length of 0 answers
// the 0 x 1 empty and a length of 1 answers the 1 x 1 one -- both of them
// before any formula runs, which is why no window here divides by N - 1
// and gets infinity.
static bool build(const std::string& name, const double& length, const Sampling& sampling,
const bool& hasParameter, const double& parameter,
ICoreMatrix& out, std::string* whyNot);
// The names build() answers, in `window`'s own order. A name outside this
// list is refused by build() with the list in the sentence.
static const std::vector<std::string>& names();
// Whether `name` takes the sampling flag, and whether it takes a
// parameter -- the evaluator asks so it can refuse `bartlett(n, 5)` with a
// sentence about bartlett rather than about the argument count.
static bool takesSampling(const std::string& name);
static bool takesParameter(const std::string& name);
// The modified Bessel function of the first kind, order zero -- kaiser's
// engine, and the one piece of new numerics this file adds. Summed from
// the power series I0(x) = sum (x/2)^(2k) / (k!)^2, whose terms are all
// POSITIVE, so there is no cancellation and the relative error stays at
// rounding for every x a window uses.
static double besselI0(const double& x);
};
};ICoreZeroPhaseFilter.h#
src/ICoreBlocks/ICoreMath/SignalProcessing/ICoreZeroPhaseFilter.h
ICoreZeroPhaseFilter#
ICoreZeroPhaseFilter.h:32 · class · 0 declaration(s)
MATLAB's filtfilt (T2.10): the same filter run forwards and then backwards, so its phase response cancels and only the squared magnitude remains.
class ICoreZeroPhaseFilter {
public:
// filtfilt(b, a, x). `x` is filtered down its COLUMNS, or along itself
// when it is a vector -- MATLAB's rule -- and the answer has x's shape.
//
// Refuses, with MATLAB's own reason, a zero leading denominator
// coefficient and a signal shorter than the pad the filter needs.
static ICoreMatrix apply(const ICoreMatrix& numerator, const ICoreMatrix& denominator,
const ICoreMatrix& x, std::string* whyNot);
// filtfilt(sos, g, x) -- second-order sections, one biquad per ROW of the
// L x 6 `sos`, with `g` the L+1 scale factors MATLAB's `tf2sos` answers
// beside them. Each stage is run through `apply`'s procedure in turn, on
// the OUTPUT of the one before it, which is what filtfilt.m's `numStage`
// loop does -- so the pad length is the stage's, not the cascade's.
//
// `g` may be empty, meaning unity: MATLAB accepts `filtfilt(sos, 1, x)`
// and a scalar gain the same way.
static ICoreMatrix applySections(const ICoreMatrix& sections, const ICoreMatrix& gains,
const ICoreMatrix& x, std::string* whyNot);
};
};