View Raw SPL
/*****************************************************************************
* *
* INVFREQZ.SPL Copyright (C) 2009 DSP Development Corporation *
* All Rights Reserved *
* *
* Author: Randy Race *
* *
* Synopsis: Computes a digital filter from a complex frequency response *
* *
* Revisions: 13 Aug 2009 RRR Creation *
* *
*****************************************************************************/
#include
#if @HELP_INVFREQZ
INVFREQZ
Purpose:
Computes digital filter coefficients from a complex
frequency response.
Syntax:
INVFREQZ(h, w, nb, na, wf, maxiter, tol)
(b, a) = INVFREQZ(h, w, nb, na, wf, maxiter, tol)
h - A series. The desired complex frequency response.
w - A series. The normalized frequency samples in
radians/s where 0 <= w <= pi and
length(w) == length(h).
nb - An integer. The order of the numerator polynomial.
na - An integer. The order of the denominator polynomial.
wf - Optional. A series, the weight at each frequency
sample. Can be an empty or a series with the same
length as h where wf > 0. Defaults to a series of
all ones (unity weighting).
maxiter - Optional. An integer, the maximum number of iterations
for the Gauss-Newton minimization option. Defaults
to 30 if specified as empty.
tol - Optional. A real, the tolerance of the gradient norm
for the iteration step. Defaults to 0.01 if specified
as empty.
Alternative Syntax:
INVFREQZ(h, w, "complex", nb, na, wf, maxiter, tol)
(b, a) = INVFREQZ(h, w, nb, na, wf, maxiter, tol)
h - A series. The desired complex frequency response.
w - A series. The normalized frequency samples in
radians/s where 0 <= w <= pi and
length(w) == length(h).
"complex" - A string. Specifies that complex filter coeffiecients
will be computed.
nb - An integer. The order of the numerator polynomial.
na - An integer. The order of the denominator polynomial.
wf - Optional. A series, the weight at each frequency
sample. Can be an empty or a series with the same
length as h where wf > 0. Defaults to a series of
all ones (unity weighting).
maxiter - Optional. An integer, the maximum number of iterations
for the Gauss-Newton minimization option. Defaults
to 30 if specified as empty.
tol - Optional. A real, the tolerance of the gradient norm
for the iteration step. Defaults to 0.01 if specified
as empty.
Returns:
A two column series where the values of column 1 are the
numerator coefficients and the values of column 2 are the
denominator coefficients.
(b, a) = INVFREQZ(h, w, nb, na, wf, maxiter, tol)
returns the numerator and denominator coefficients as two
separate series.
If "complex" is specified, complex coefficients are returned.
Example:
b = {1, 2, 3, 2, 3}
a = {1, 2, 3, 2, 1, 4}
(h, w) = freqz(b, a, 64);
(b1, a1) = invfreqz(h, w, 4, 5)
b1 == {1, 2, 3, 2, 3}
a1 == {1, 2, 3, 2, 1, 4}
The variables h and w contain the complex frequency
response and frequency values of the system to model.
Although INVFREQZ computes the filter coefficients that
match the original values, because the denominator contains
poles outside the unit circle, the system is unstable.
Example:
(b2, a2) = invfreqz(h, w, 4, 5, {}, 20)
b2 == {0.249368, 0.301802, -0.001658, 0.063888, 0.121041}
a2 == {1.000000, -0.794051, 0.575584, 0.903462, -0.697508, 0.479046}
By specifying an iteration value, INVFREQZ computes the
filter coefficients using the Gauss-Newton minimization
method. The computed system is stable (all denominator
poles are on or inside the unit circle).
Remarks:
INVFREQZ computes digital filter coefficients given a
complex frequency response and normalized frequency values
for the following discrete system:
B(z) b[1] + b[2]*z^(-1) + ... + b[n+1]*z^(-n)
H(z) = ---- = ----------------------------------------
A(z) a[1] + a[2]*z^(-1) + ... + a[m+1]*z^(-m)
The parameters nb and na determine the order of the numerator
and denominator polynomials such that nb+1 and na+1 coefficients
will be returned.
If maxiter is not specified, the coefficients are computed
by minimizing:
sum(b - h*a)^2 * wf
If maxiter is specified as empty or an integer, the coefficients
are computed using the Gauss-Newton method to minimize:
sum(b/a - h)^2 * wf
The input frequencies, w, are specified in radians/seconds
such that
w = 2*pi*f
where f is frequency in Hertz.
If "complex" is not specified, real coefficients are
computed. In this case, the normalized frequency values
must range from 0 to pi and the response of the filter at
negative frequencies is automatically set to conj(h) to
maintain proper symmetry.
If "complex" is specified, complex filter coefficients are
computed and the normalized frequency samples can range
from -pi to pi and no frequency domain symmetry is enforced.
The first coefficient of a, the denominator polynomial, is 1.0.
See Also:
Filteq
Freqs
Freqz
Gimpulse
Grpdelay
Impz
Invfreqs
Residue
Residuez
Sfreq
Zfreq
#endif
/* compute digital filter coefficients from complex frequency response */
invfreqz(h, w, complex, nb, na, wf, maxiter, tol, pf)
{
local a, b, iter;
/* parse input args */
(h, w, complex, nb, na, wf, maxiter, tol, pf, iter) = invfreqz_parse_args(h, w, complex, nb, na, wf, maxiter, tol, pf);
/* the computional routine */
(b, a) = invfreq(h, w, complex, nb, na, wf, maxiter, tol, pf, 'z', iter);
if (outargc > 1)
{
return(b, a);
}
else
{
/* combined series */
return(ravel(b, a));
}
}
invfreqz_parse_args(h, w, complex, nb, na, wf, maxiter, tol, pf)
{
local iter;
switch (argc)
{
case 0:
error("invfreqz - input gain a frequencies required");
break;
case 1:
w = xvals(h);
case 2:
complex = "";
case 3:
nb = {};
case 4:
na = {};
case 5:
wf = {};
case 6:
maxiter = {};
case 7:
tol = {};
case 8:
pf = {};
case 9:
iter = {};
case 10:
break;
default:
error("too many arguments for invfreqz");
break;
}
/* check complex flag */
if (not(isstring(complex)))
{
pf = tol;
tol = maxiter;
maxiter = wf;
wf = na;
na = nb;
nb = complex;
complex = "";
iter = argc > 5;
}
else
{
iter = argc > 6;
}
if (isempty(nb))
{
nb = 1;
}
if (isempty(na))
{
na = 1;
}
if (isempty(wf))
{
wf = ones(length(w), 1);
}
return(h, w, complex, nb, na, wf, maxiter, tol, pf, iter);
}