View Raw SPL
/*****************************************************************************
* *
* RESAMPLE.SPL Copyright (C) 2008-2018 DSP Development Corporation *
* All Rights Reserved *
* *
* Author: Randy Race *
* *
* Synopsis: Resamples a series to an arbitrary sample rate *
* *
* Revisions: 25 Jul 2008 RRR Creation *
* 20 Jan 2012 RRR rate == 0 implies no Y interpolation *
* 14 Feb 2017 RRR nearest neighbor offset fix *
* 23 May 2018 RRR sinxinterp and sinxinterpxy *
* 4 Nov 2018 RRR sincinterp and sincinterpxy *
* *
*****************************************************************************/
#include
#if @HELP_RESAMPLE
RESAMPLE
Purpose: Resamples a series to an arbitrary sample rate.
Syntax: RESAMPLE(series, rate, method)
series - A series. The input series to resample.
rate - A real, the new sample rate.
method - Optional. A string. The interpolation method.
"linear" : use linear interpolation (default)
"spline" : use cubic spline
"sinx" : use periodic sinx/sin
"none" : return nearest neighbor (no interpolation)
Alternate Syntax:
RESAMPLE(xseries, yseries, rate, method)
xseries - A series. The input X series to resample.
yseries - A series. The input Y series to resample.
rate - A real, the new sample rate.
method - Optional. A string. The interpolation method.
"linear" : linear interpolation (default)
"spline" : cubic spline
"sinx" : periodic sinx/sin
"sinc" : non-periodic sin(x)/x
"none" : nearest neighbor (no interpolation)
Returns: An interval series where the Y values are sampled at the given
rate.
Example:
W1: {1, 2, 3, 2}
W2: resample(w1, 3)
W3: W1;overp(w2,lred);setsym(14,2);setplotstyle(1,2)
W1 contains the source series with a sample rate of 1.
W2 linearly resamples W1 to a rate of 3.0 and produces
a 10 point series. The time duration of W2 is the same as
W1.
W3 displays the original series with the resampled values
overplotted as red circles.
Example:
W4: resample(w1, 3, "spline")
Same as above except cubic spline interpolation is performed.
Example:
W1: {1, 1, 2, 2, 3}
W2: {1, 2, 3, 4, 5}
W3: xy(w1, w2)
W4: resample(w3, 0)
W3 contains an XY series with duplicate X values. W4 creates
an interval series by preserving the original Y values
and spanning the X values from the minimum X value to the
maximum X value. No interpolation is performed.
Remarks:
If RATE is set to 0.0, the original Y values are preserved
and the resulting sample rate is based on the minimum and
maximum X values.
See Also:
Interp
Spline
Xyinterp
Xylookup
#endif
/* resample to a given sample rate */
ITERATE resample(x, y, rate, method)
{
local intv, yout, dx, xi, hunits;
/* parse args */
if (argc < 4)
{
if (not(isarray(x))) resample_error_mes1();
if (argc < 3)
{
if (argc < 2)
{
/* resample(y) */
y = refseries(x);
rate = -1;
method = "linear";
intv = 1;
}
else
{
if (isarray(y))
{
/* resample(x, y) */
rate = -1;
method = "linear";
intv = 0;
}
else
{
if (isstring(y))
{
/* resample(y, method) */
method = y;
y = refseries(x);
rate = 1.0 / castreal(xyinterp(y, -1, 1));
intv = 1;
}
else
{
/* resample(y, rate) */
rate = castreal(y);
y = refseries(x);
method = "linear";
intv = 1;
}
}
}
}
else
{
/* resample(x, y, rate) or resample(y, rate, method) */
if (isarray(y))
{
if (isstring(rate))
{
/* resample(x, y, method) */
method = rate;
rate = 1.0 / castreal(xyinterp(x, y, -1, 1));
intv = 0;
}
else
{
/* resample(x, y, rate) */
rate = castreal(rate);
method = "linear";
intv = 0;
}
}
else
{
/* resample(y, rate, method) */
if (isstring(y))
{
method = y;
}
else
{
method = caststring(rate);
rate = castreal(y);
}
y = refseries(x);
intv = 1;
}
}
}
else
{
if (not(isarray(x)) && not(isarray(y)))
{
resample_error_mes1();
}
if (isstring(rate))
{
/* swap */
(rate, method) = (method, rate);
}
rate = castreal(rate);
method = caststring(method);
intv = 0;
}
if ((dx = rate) != 0.0)
{
dx = 1 / rate;
}
else
{
/*
* rate == 0, no interpolation, use X range,
* xyinterp will do the job
*/
method = "linear";
}
/* now do calculation */
switch (tolower(method))
{
case "linear":
if (iscomplex(y))
{
yout = intv ? resample_xyinterp_complex(y, dx) : resample_xyinterp_complex(x, y, dx);
}
else
{
yout = intv ? xyinterp(y, dx) : xyinterp(x, y, dx);
}
break;
case "sinx":
case "periodic":
case "fft":
yout = intv ? sinxinterp(y, rate) : sinxinterpxy(x, y, rate);
break;
case "sinc":
yout = intv ? sincinterp(y, rate) : sincinterpxy(x, y, rate);
break;
case "spline":
yout = intv ? csinterp(y, rate) : csinterpxy(x, y, rate);
break;
case "none":
/* nearest neighbor */
if (intv)
{
if (isxy(y))
{
x = xvals(y);
xi = min(x)..dx..max(x);
}
else
{
xi = xoffset(y)..dx..((length(y) - 1) * deltax(y) + xoffset(y));
}
yout = xylookup(y, xi, "none");
hunits = gethunits(y);
}
else
{
xi = min(x)..dx..max(x);
yout = xylookup(x, y, xi, "none");
hunits = getvunits(x);
}
/* convert to interval series */
yout = yvals(yout);
setxoffset(yout, min(xi));
setdeltax(yout, dx);
sethunits(yout, hunits);
break;
default:
error(sprintf('%s - unknown method "%s"', __FUNC__, method));
break;
}
return(yout);
}
resample_error_mes1()
{
error("resample - input series required");
}
/* complex xyinterp */
resample_xyinterp_complex(x, y, dx)
{
local yr, yi, yo;
if (argc == 3)
{
yr = xyinterp(x, real(y), dx);
yi = xyinterp(x, imag(y), dx);
}
else
{
yr = xyinterp(real(x), y);
yi = xyinterp(imag(x), y);
}
yo = yr + i * yi;
return(yo);
}