View Raw SPL
/*****************************************************************************
* *
* SINFIT3.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_SINFIT3
SINFIT3
Purpose: Fits y(x) = A * cos(2*pi*F*x + B) + C using least squares
Syntax: SINFIT3(s, freq, mode)
s - A series, the input.
freq - Optional. A real, the known frequency in hertz.
Defaults to -1, estimate frequency.
mode - specifies output in one of the two following forms:
0 - returns the fitted curve and optionally the coefficients
and rms error for the equation:
y(t) = A * cos(2 * pi * Freq * t) +
B * sin(2 * pi * Freq * t) + C
If the output is directed to variables:
(fit, coef) = sinfit3(s, freq, mode)
Mode 0 returns coefficients in the form
coef = {A, Freq, B, C}
1 - (default) returns the fitted curve and optionally the
coefficients and rms error for the equation:
y(t) = A * cos(2 * pi * Freq * t + B) + C
If the output is directed to variables:
(fit, coef) = sinfit3(s, freq, mode)
Mode 1 returns coefficients in the form
coef = {A, Freq, B, C}
Returns: A series, the fitted sine curve.
(fit, coef) = sinfit(s)
returns both the fit and the coefficients as a series.
(fit, coef, rmserr) = sinfit(s)
returns the fit, coefficients and RMS error of the fit.
Example:
W1: 5 + 3 * gsin(1000, .001, 4, pi) + gnorm(1000, .001)
W2: sinfit3(w1);overp(w1, lred)
Overplots the original data with the calculated sin fit.
Example:
(fit, coef) = sinfit3(w1)
fit is the same series as in W2
coef == {3.0012, 3.9999, 1.5708, 5.000}
Example:
(fit, coef, rmserr) = sinfit3(w1, 4)
fit contains the fitted data
coef == {3.0012, 4.0000, 1.5706, 5.000}
rmserr == 0.0999
Remarks:
SINFIT3 uses the least squares method to fit a sinusoid to
a known frequency as per the IEEE 1241 Standard. If the
input frequency is unspecified, the frequency is estimated
by locating the dominant frequency of the FFT and using
analytical rectangular window interpolation to determine
the precise frequency. The FFT interpolation handles
estimated frequencies that do not lie on a unique FFT bin
frequency.
The frequency coefficient, F = coef[2] is in Hertz and the
phase term, C = coef[3] is in radians.
SINFIT3 uses LSINFIT to perform the least squares fit as per
the IEEE 1241 Standard.
See SINFIT4 for an implementation of the 4 term iterative
sinusoid fit as outlined in the IEEE 1241 Standard.
See SINFIT for a similar routine that uses Hanning window
interpolation and returns the frequency coefficent in radians.
See Also:
FFT
Linfit
Lsinfit
Sinfit
Sinfit4
Sintrend
Trend
References:
IEEE Std 1241-2010 Annex B.1 "An algorithm for three-parameter
(known frequency) least-squares fit to sine-wave data"
#endif
/* least squares fit y = A * cos(2*pi*F*x + B) + C */
ITERATE sinfit3(s, freq, mode)
{
local f, n, Fc, fmag, coef, coef3, n0, rmserr;
local U, V, U1, V1, Kopt, Z1, Z2, lambda, fser;
if (argc < 3)
{
if (argc < 2)
{
if (argc < 1) error("sinfit - input series required");
/* estimate frequency */
freq = -1;
}
/* cosine form, A * cos(2*pi*F*x + B) + C */
mode = 1;
}
/* degenerate case */
if (length(s) <= 1)
{
if (freq < 0)
{
freq = 0;
}
return({max(s)}, {0, max(s), freq, 0});
}
/* check if estimate frequency */
if (freq <= 0)
{
if (isxy(s))
{
/* interpolate XY to interval series */
ser = xyinterp(s);
}
else
{
ser = refseries(s);
}
/* series length */
n = length(ser);
/* non-windowed FFT */
f = fft(demean(ser), n);
/* first half */
Fc = extract(f, 1, int(n/2));
/* normalized magnitude */
fmag = mag(Fc) / length(ser);
/* index of primary frequency */
w = maxidx(fmag);
/* use left or right side of max for interpolation */
if (w > 1)
{
if (fmag[w - 1] > fmag[w + 1])
{
w--;
}
w--;
}
/* perform rectangular window interpolation */
n0 = 2 * pi / n;
U = real(Fc[w + 1]);
V = imag(Fc[w + 1]);
U1 = real(Fc[w + 2]);
V1 = imag(Fc[w + 2]);
Kopt = (sin(n0 * w) * (V1 - V) + cos(n0 * w) * (U1 - U)) / (U1 - U);
Z1 = V * (Kopt - cos(n0 * w)) / sin(n0 * w) + U;
Z2 = V1 * (Kopt - cos(n0 * (w + 1))) / sin(n0 * (w + 1)) + U1;
lambda = acos((Z2 * cos(n0 * (w + 1)) - Z1 * cos(n0 * w)) / (Z2 - Z1)) / n0;
lambda = acos((Z2 * cos(n0 * (w + 1)) - Z1 * cos(n0 * w)) / (Z2 - Z1));
lambda = lambda / n0;
/* freq in hertz */
freq = lambda * deltax(fmag);
}
/* three term least squares fit */
(fit, coef3, rmserr) = lsinfit(s, freq, mode);
if (outargc > 1)
{
/* arrange coefficients */
coef = zeros(4, 1);
/*
* mode == 0, A * cos(2*pi*F*x) + B * sin(2*pi*F*x) + C
*
* mode == 1, A * cos(2*pi*F*x + B) + C
*/
coef = {coef3[1], freq, coef3[2], coef3[3]};
if (outargc > 2)
{
return(fit, coef, rmserr);
}
else
{
return(fit, coef);
}
}
else
{
return(fit);
}
}