nutation_dist.m

Nutation frequency distribution from a nutation curve measured with the same coil used for excitation and detection. The trace is normalised to unit maximum modulus, and the complex noise variance is estimated from the median of the squared second difference of the samples, capped at 5% of the peak power and floored at the machine-precision level. The receiver phase line is obtained from the argument of the sum of the squares of the samples, an estimate that carries a pi ambiguity; both signs are tried. A raised-cosine fade-out window and an eight-fold zero-filled FFT give a sine-transform power spectrum, in which the noise-limited support is the set of bins above ten noise bins, widened by two frequency steps on either side, doubled about its centre, and clipped to the Nyquist interval. The reconstruction grid holds between 160 and 400 points at approximately twice the Fourier resolution. The fitting kernel is sin((t+t0)*freq) multiplied by the dimensionless reciprocity reception weight freq/freq_hi, and the distribution is obtained by non-negative least squares on the data stacked with a sqrt(lambda)-scaled second-difference penalty, with the receiver phase re-estimated by four alternation cycles on each of the two phase branches. The time origin error t0 is selected from nine trial shifts spanning plus or minus two sampling intervals, and then refined on the surrounding interval with five more. The fitted masses are clipped at zero and divided by their trapezoidal integral.

Syntax

    [freq,distr]=nutation_dist(curve,dt,lambda)

Parameters

    curve  - nutation curve, a row or column vector; either
             the complex X+iY output of a quadrature receiver,
             or a real phase-corrected trace

    dt     - sampling interval in seconds

    lambda - second-derivative Tikhonov regularisation para-
             meter, a non-negative real scalar; the curve is
             normalised to unit maximum modulus and the fit-
             ting kernel is dimensionless, so lambda is a di-
             mensionless number of the order of the ratio of
             the squared Frobenius norms of the kernel and of
             the second difference matrix; zero switches the
             regularisation off

Outputs

    freq   - nutation frequency grid in rad/s, a column vector

    distr  - non-negative nutation frequency density in inverse
             rad/s, normalised to unit integral over freq

Notes

When one coil both excites and detects, the reciprocity principle makes the detected amplitude of every isochromat proportional to its own nutation frequency, and only the sine component of the nutation is observable. The curve is therefore modelled as an unknown complex receiver scale times

    s(t)=integral(distr(freq)*freq*sin(freq*(t+t0)),d freq)

and the reception weight is divided out, so that the returned density is the true nutation frequency distribution. A real receiver scale is a valid special case, and a phase-corrected real trace is therefore accepted. The frequency support, receiver phase, and sub-sample time shift are selected from the supplied trace; the returned density is therefore a stable estimate rather than a raw finite-time transform. The reconstruction runs on twice the noise-limited support width, centred on the same band and clipped to the Nyquist interval, so that the margins of the distribution are visible rather than truncated.

See also

tikhonov.m, tikhoind.m, apodisation.m, Kernel utilities

Version 2.13, authors: Ilya Kuprov