View Raw SPL
/*****************************************************************************
* *
* SINFIT4.SPL Copyright (C) 2012 DSP Development Corporation *
* All Rights Reserved *
* *
* Author: Randy Race *
* *
* Synopsis: Fits y(x) = A * cos(2*pi*F*x + B) + C using least squares *
* *
* Revisions: 15 May 2012 RRR Creation *
* *
*****************************************************************************/
#if @HELP_SINFIT4
SINFIT4
Purpose: Fits y(x) = A * cos(2*pi*F*x + B) + C by iterative least squares.
Syntax: SINFIT4(s, freq, niter, tol)
s - A series, the input.
freq - Optional. A real, the known frequency in hertz.
Defaults to -1, estimate frequency.
niter - Optional. An integer, the maximum number of least squares
iterations if no convergence. Defaults to 30.
tol - Optional. A real, the relative convergence tolerance.
Defaults to 1e-15.
Returns: A series, the fitted sine curve.
(fit, coef) = sinfit4(s)
returns both the fit and the coefficients as a series.
(fit, coef, rmserr) = sinfit4(s)
returns the fit, coefficients and RMS error of the fit.
Example:
W1: 5 + 3 * gsin(1000, .001, 4, pi) + gnorm(1000, .001)
W2: sinfit4(w1);overp(w1, lred)
Overplots the original data with the calculated sin fit.
Example:
(fit, coef) = sinfit4(w1)
fit is the same series as in W2
coef == {2.996476, 4.009087, 1.549703, 4.999965}
Example:
(fit, coef, rmserr) = sinfit4(w1, 4)
fit contains the fitted data
coef == {2.006563, 4.008975, 1.5749337, 4.999968}
rmserr == 0.0999
Remarks:
SINFIT4 uses the iterative least squares method to fit a
sinusoid with a starting frequency as per the IEEE 1241
Standard such that:
y(t) = A * cos(2*pi*F*t + theta) + C
If the starting frequency is unspecified, the frequency is
estimated by using SINFIT3. SINFIT3 is also used as the
starting point for the iterative least squares fit. Given
the specified or estimated starting frequency, the
frequency is adjusted to minimize the least squares error
to the tolerance value.
The frequency coefficient, F = coef[2] is in Hertz and the
phase term, C = coef[3] is in radians.
The starting frequency, if specified, may differ from the
resulting fitted frequency. See SINFIT3 to perform a fit
that preserves a known frequency. SINFIT3 is an
implementation of the 3 term sinusoid fit as outlined in
the IEEE 1241 Standard.
See Also:
FFT
Linfit
Lsinfit
Sinfit3
Sintrend
Trend
References:
IEEE Std 1241-2010 Annex B.2 "An algorithm for four-parameter
least-squares fit to sine-wave data"
#endif
/* fit y = A * cos(2*pi*F*t + B) + C */
ITERATE sinfit4(s, freq, niter, tol)
{
local M, fit, coef, rmserr, t, werror, w, k, oiter;
local A, F, B, C;
if (argc < 4)
{
if (argc < 3)
{
if (argc < 2)
{
if (argc < 1) error("sinfit4 - input series required");
/* estimate frequency */
freq = -1;
}
/* default iterations */
niter = -1;
}
/* default tolerance */
tol = -1;
}
/* degenerate case */
if (length(s) <= 1)
{
if (freq < 0)
{
freq = 0;
}
return({max(s)}, {0, max(s), freq, 0});
}
if (niter < 0)
{
/* default iterations */
niter = 30;
}
if (tol < 0)
{
/* default tolerance */
tol = 1e-15;
}
loop (k = 1..niter)
{
(coef, oiter, t) = sinfit4_iter(s, freq, niter, tol);
if (oiter < niter)
{
break;
}
seedrand(0);
/* increase tolerance */
tol *= 2;
/* adjust frequency */
freq = freq + freq / 80 * randn();
}
/* coefficients, A * cos(2*pi*F*x + B) + C */
w = coef[4];
A = sqrt(coef[1]^2 + coef[2]^2);
F = w / (2 * pi);
B = atan2(-coef[2], coef[1]);
C = coef[3];
/* fit */
fit = A * cos(w * t + B) + C;
if (outargc > 1)
{
/* coefficients */
coef = {A, F, B, C};
if (outargc > 2)
{
/* rms error */
rmserr = sqrt(mean((s - fit)^2));
return(fit, coef, rmserr);
}
else
{
return(fit, coef);
}
}
else
{
return(fit);
}
}
sinfit4_iter(s, freq, niter, tol)
{
local fit, t, coef, ocoef, r, A, B, M, k, wt;
/* estimate initial parameters via 3 term least squares fit */
(fit, coef) = sinfit3(s, freq, 0);
/* amplitudes from 3 parameter fit */
A = coef[1];
B = coef[3];
if (isxy(s))
{
r = rate(fit);
s = yvals(s);
}
else
{
r = rate(s);
}
/* time values */
t = xvals(fit);
/* angular frequency */
w = 2 * pi * coef[2];
/* approximation tolerance in radian frequency */
werror = 2 * pi * tol * r;
coef = {A, B, coef[2], 0.0};
/* least squares iteration on w */
loop (k = 1..niter)
{
/* adjust radian frequency */
w += coef[4];
wt = w * t;
ocoef = coef;
/* least squares design matrix */
M = ravel(cos(wt), sin(wt), ones(length(t), 1), -A * t * sin(wt) + B * t * cos(wt));
/* solve */
coef = linfit(M, s);
/* check tolerance */
if (abs(coef[4]) < werror)
{
break;
}
/* new amplitudes */
A = coef[1];
B = coef[2];
}
/* new frequency */
coef[4] = w;
return(coef, k, t);
}