View Raw SPL
/*****************************************************************************
* *
* FILTEQ.SPL Copyright (C) 1999-2015 Development Corporation *
* All Rights Reserved *
* *
* Author: Randy Race *
* *
* Synopsis: Evaluates a Linear Constant Coeffcient Difference Equation *
* *
* Revisions: 22 Mar 1999 RRR Creation *
* 1 Aug 2001 RRR Optional Z transform form *
* 10 Oct 2003 RRR Added optional xi initial conditions *
* 18 Nov 2009 RRR Cast inputs to series *
* 15 Mar 2011 RRR FIR optimization and final conditions *
* 5 Jun 2015 RRR Faster DIFFEQ implementation *
* *
*****************************************************************************/
#if @HELP_FILTEQ
FILTEQ
Purpose: Evaluates a Linear Constant Coefficient Difference Equation
Syntax: FILTEQ(b, a, x, yi, xi)
b - a series, x[n] (FIR or zero) coefficients
a - a series, y[n] (IIR or pole) coefficients,
if a[1] == 1, the coefficients are assumed to be
in standard Z transform form
x - a series, the input
yi - an optional series, Y initial conditions
xi - an optional series, X initial conditions
Returns: A series.
(y, iirzf, firzf) = FILTEQ(b, a, x) returns the filtered
series and final conditions.
Example:
x = {1, 0, 0, 0, 0}
y = filteq({1, -0.5}, {0.8, -0.2, -0.5}, x)
y == {1.0, 0.3, 0.04, -0.528, -0.5804}
y is the output produced by the difference equation:
y[n] = x[n] - 0.5*x[n-1] + 0.8*y[n-1] - 0.2*y[n-2] - 0.5*y[n-3]
Example:
W1: impulse(1, 1024)
W2: filteq({1, -0.5}, {0.8, -0.2, -0.5}, W1)
W3: spectrum(w2)
calculates 1024 samples of the impulse response of the
above difference equation. The Spectrum in W3 shows a
resonate peak at approximately 0.127 Hertz.
Example:
x = {1, 0, 0, 0, 0}
y = filteq({1, -0.5}, {1, -0.2, -0.5}, x)
y == {1.0, -0.3, 0.44, -0.062, 0.2076}
Because a[1] == 1, y is the output produced by the difference
equation:
y[n] = x[n] - 0.5*x[n-1] + 0.2*y[n-1] + 0.5*y[n-2]
or the equivalent Z transform:
-1
Y(z) 1 - 0.5*z
H(z) = ---- = -----------------------
-1 -2
X(z) 1 - 0.2*z - 0.5*z
Example:
b = {1, -0.5};
a = {1, 0.5, -0.5};
W1: gnorm(200, 1);
W2: filteq(b, a, w1);
W3: extract(w1, 1, int(length(w1)/2));
W4: extract(w1, 1 + int(length(w1)/2), -1);
W5: (y, iirzf, firzf) = filteq(b, a, w3);y;
W6: filteq(b, a, w4, iirzf, firzf);
W7: concat(w5, w6);
W8: w2 - w7;
W1 contains 200 samples of random noise. W2 filters the series
with a finite difference equation. W3 and W4 divide the series
in W1 into two segments. W5 filters the first segment and returns
the final conditions. W6 filters the second segment setting the
initial conditions to the final conditions of W5. W7 combines
the output results and W8 shows the two processes return the same
filtered output series.
Remarks:
FILTEQ evaluates a linear constant coefficient difference
equation of the following form:
y[n] = a1*y[n-1] + a2*y[n-2] + a3*y[n-3] + ... + aN*y[n-N] +
b0*x[n] + b1*x[n-1] + b2*x[n-2] + ... + bM*x[n-M]
for n >=1
Initial conditions, y[0], y[-1], ... y[-L] are specified by
the optional series yi. If not specified, yi is assumed to
be 0.
Initial conditions, x[0], x[-1], ... x[-K] are specified by
the optional series xi. If not specified, xi is assumed to
be 0.
If a[1] == 1, the coefficients are assumed to be in standard
Z transform form:
-1 -Q
Y(z) b[1] + b[2]*z + ... + b[Q+1]*z
H(z) = ---- = -------------------------------------
-1 -P
X(z) 1 + a[2]*z + ... + a[P+1]*z
See Also:
Fir
Iir
Zfreq
#endif
/* LCC difference equation */
ITERATE filteq(b, a, x, yi, xi)
{
local y, iirzf, firzf;
/* check inputs */
if (argc < 5)
{
if (argc < 4)
{
if (argc < 3) error("filteq - b, a, and input series required");
yi = {};
}
xi = {};
}
/* cast to series */
b = {b};
a = {a};
x = {x};
/* check for empty series */
if (length(b) < 1)
{
b = {0};
}
if (length(a) < 1)
{
a = {0};
}
else if (filteq_isequal(a[1], 1.0))
{
/* Z transfer form - convert to difference equation */
if (length(a) == 1)
{
a = {0};
}
else
{
a = -extract(a, 2, -1);
}
}
/* difference equation - handles empty initial conditions */
if (any(a))
{
/* both iir and fir terms */
if (outargc > 0)
{
if (length(xi) > 0 || outargc > 2)
{
/* return final conditions - old method */
y = iir(fir(x, b, {xi}), a, {yi});
/* iir final conditions */
iirzf = rev(extract(y, length(y) - length(a) + 1, -1));
/* fir final conditions */
firzf = rev(extract(x, length(y) - length(b) + 1, -1));
return(y, iirzf, firzf);
}
else
{
/* faster diffeq */
if (outargc > 1)
{
firzf = {};
(y, iirzf) = diffeq(b, {1, -a}, x, {yi});
return(y, iirzf, firzf);
}
else
{
return(diffeq(b, {1, -a}, x, {yi}));
}
}
}
else
{
if (length(xi) > 0)
{
/* old method */
return(iir(fir(x, b, {xi}), a, {yi}));
}
else
{
return(diffeq(b, {1, -a}, x, {yi}));
}
}
}
else
{
/* purely fir */
if (outargc > 0)
{
/* return final conditions */
y = fir(x, b, {xi});
/* iir final conditions */
iirzf = {};
/* fir final conditions */
firzf = rev(extract(x, length(y) - length(b) + 1, -1));
return(y, iirzf, firzf);
}
else
{
return(fir(x, b, {xi}));
}
}
}
/* are two values equal within a tolerence */
filteq_isequal(d1, d2, tol)
{
local eq;
if (argc < 3)
{
tol = 1e-8;
}
eq = abs(d1 - d2) < abs(d1 * tol);
return(eq);
}