View Raw SPL
/*****************************************************************************
* *
* IMPS.SPL Copyright (C) 2015-2016 DSP Development Corporation *
* All Rights Reserved *
* *
* Author: Randy Race *
* *
* Synopsis: Impulse response of an S transform *
* *
* Revisions: 5 Aug 2015 RRR Creation *
* 1 Jun 2016 RRR repeated pole support *
* *
*****************************************************************************/
#if @HELP_IMPS
IMPS
Purpose: Calculates the impulse response of an S transform.
Syntax: IMPS(b, a, t)
(h, t) = IMPS(b, a, t)
b - A series, the numerator coefficients in descending
powers of s.
a - A series, the denominator coefficients in descending
powers of s.
t - Optional. A series or integer, the time values. If an
integer or not specified, the time duration is computed
based on the exponential time constant and t values are
returned where t defaults to 512. If a series, the t
values are used directly.
Returns: An XY series, the impulse response of the system.
(h, t) = IMPS(b, a, t) returns the impulse response and time
values as two separate series.
Example:
1
Given H(s) = -------
s + 0.5
W1: imps({1}, {1, 0.5})
returns 512 samples of the series h(t) = exp(-0.5 * t) in W1.
Example:
1 1
Given H(s) = ------------ = --------------
s^2 + 4s + 3 (s + 3)(s + 1)
H(s) has the following partial fraction expansion:
-0.5 0.5
H(s) = ----- + -----
s + 3 s + 1
The impulse response for t >= 0 can be found by inspection:
h(t) = -0.5 * exp(-3 * t) + 0.5 * exp(-t)
W1: imps({1}, {1, 4, 3})
W2: -0.5 * exp(-3*xvals(w1)) + 0.5 * exp(-xvals(w1))
W3: w1 - w2
W1 displays 512 samples of the calculated response.
W2 contains 512 samples of the analytical impulse response
using the same time values as W1.
W3 confirms that the difference is negligible.
Remarks:
The input series represent the terms of the rational polynomial
H(s) = b(s)/a(s) where:
N = length(b) and M = length(a):
b(s) b[1] s^(N-1) + b[2] s^(N-2) + ... + b[N]
H(s) = ------ = ----------------------------------------
a(s) a[1] s^(M-1) + b[2] s^(M-2) + ... + a[M]
The impulse response is calculated from the partial fraction
expansion of the system equation.
If there are no output arguments, the result is displayed in
the current window.
See Also:
Impz
Residue
Sfreq
Splane
#endif
/* continuous filter impulse response */
imps(b, a, t)
{
local r, p, k, f, j, dur, N, s, m, n;
if (argc < 3)
{
if (argc < 2)
{
if (argc < 1) error("imps - input series required");
a = {};
}
t = {};
}
/* find residues and poles */
(r, p, k) = residue(b, a);
if (length(k) > 1)
{
error("imps - number of poles must be >= number of zeros");
}
if (length(t) < 2)
{
if (isscalar(t))
{
N = (t > 0) ? int(t) : 512;
}
else
{
N = 512;
}
/* estimate time duration */
dur = imps_duration(r, p, N);
t = linspace(0.0, dur, N);
/* units */
setvunits(t, "s");
}
else
{
if (getvunits(t) == "No Units")
{
setvunits(t, "s");
}
}
n = 0;
m = 1;
f = (length(k) > 0) ? k : 0;
/* impulse response of each single pole section */
loop (j = 1..length(r))
{
s = r[j] * exp(p[j] * t);
if (j > 1)
{
if (p[j] == p[j-1])
{
/* repeated pole */
n++;
m *= n;
s *= t^n / m;
}
else
{
n = 0;
m = 1;
}
}
f += s;
}
/* check if real */
if ((isreal(b) && isreal(a)) || all(abs(imag(f)) < eps("float")))
{
f = real(f);
}
if (outargc > 1)
{
return(f, t);
}
else
{
return(xy(t, f));
}
}
/* estimated impulse response duration duration */
imps_duration(rr, pp, N)
{
local tc, dur, tset = 0.005, dt, p;
/* remove zeros */
p = pp[find(not(pp == 0))];
/* exponential time constant */
tc = log(tset) / real(p);
/* duration */
dur = max(abs(tc));
/* sample rate */
dt = min(min(tc) / 5, max(tc) / 50);
dt = min(dt, pi / (10 * max(mag(pp))));
/* constrain to number of samples */
if (dur / dt > N)
{
dur = dt * N;
}
return(dur);
}