In article <o4una7l1d9mser16dp1el1vp4ernmi5726@4ax.com>,
gjberchin@chatter.net.invalid says...
>Thanks, Robert. Very interesting. I have some questions and comments:
>
>>First we frequency-warp the desired z-plane
>>magnitude function so that it can be optimized as a quotient of two
>>polynomials with no trig terms. Basically, the frequency warp uses an
>>"inverse bilinear transformation." In other words, this creates an s-plane
>>magnitude response that will produce the desired z-plane function once
>>the normal bilinear transformation is applied to the s-plane poles and
>>zeros that we find in the optimization.
>
>So you perform this initial "unwarping" with the knowledge that your final
step
>will be warped by a bilinear transform. Have you experimented with other s-
>z
>mapping methods and their attendant warping (if any)? For example, if you
know
>that your final step will be a "matched pole-zero" mapping, then IIRC there
will
>be no warping at all.
I have not experimented with other mappings, although I suppose that you
could use any of the common ones like impulse-invariant, step-invariant, or
matched-z (as you said). Once nice thing about the inverse bilinear warp is
that if you uniformly frequency-sample the frequency axis after the warp, it
automatically concentrates more grid points close to what will be Nyquist
after the final "unwarp," which is where the optimizer tends to have the most
trouble. Another nice thing is that if the inner optimization has a minimax
error curve, the minimax property will be preserved in the final result
because of the properties of the bilinear transform -- it doesn't change the
magnitudes; it just moves them to different frequencies.
The other main goal of the warp is to allow the optimizer to operate on a
ratio of polynomials containing no trig or exponential functions, which
allows one to choose from a wide range of well established optimization
techniques without gumming up the mathematics. In addition, I suspect that it
speeds up the optimizer because it doesn't need to calculate trig functions.
One thing that I didn't mention in my post is that although in principle, the
number of z-plane poles and zeros can be different in the optimization, the
technique is mainly useful if the number of z-plane poles and zeros is forced
to be the same because this allows the magnitude response at Nyquist to be
something other than zero or infinite after the final bilinear unwarp. This
makes life much easier on the optimizer because it prevents a very large
magnitude change from happening in the frequency range between Nyquist and
upper frequency range of the optimization.
>
>>Next, we make a second warp of the frequency axis by changing the
>>frequency variable from w^2 to u,
>
>To state an implicit step explicitly, I think you're saying that you compute
the
>magnitude-squared expression and then apply the change of variables.
Exactly correct.
>
>Again, I wonder: if you used a non-warping s->z mapping method, AND you knew
the
>exact Laplace domain transfer function for the prototype frequency response,
>could you simply use the magnitude-squared form of that transfer function at
>this point?
The main point of the change of variables is to minimize the degree of the
polynomials in the optimization. Rootfinding occurs as part of the Remez
update (in fact, the rootfinding has to be applied to a mildly nonlinear
equation so there is an "inner loop" in the Remez update that deals with
this; see the Acton reference for more detail), so I thought that minimizing
the degree of the polynomials would help the rootfinder. There are some
"interesting" numerical stability issues that arise in the Remez update, so
it is important to choose a good rootfinder (I used SVD) and to give it all
of the help one can.
>
><snip minimax implementation details>
>
>>The rest of the algorithm just requires taking the resulting ratio of
>>polynomials and frequency-mapping the poles and zeros inversely to
>>the way the original mapping was done. This is done in two steps: the
>>first replaces u by s^2, and the second applies the normal bilinear
>>transform to the poles and zeros in s. This gets us to the z plane and the
>>final poles and zeroes of our optimization.
>
>My suspicion is that, if you used a non-warping mapping method, you could
apply
>your technique directly to a *known* Laplace domain transfer function. Thus
you
>would have a direct s->z magnitude function approximation.
I would have to think this through; you may be correct. Because it is based
on a bilinear warp, my technique preserves the minimax error in the
optimization through the two steps referenced in my quote. So if you want a
true minimax error curve, my method is a surefire way to obtain it. I'm not
sure if this would be true of a non-warping method; again, it's something I
would have to think through.
>
>(Repeating the example cited above, your method would produce different
results
>than simple "matched pole-zero" mapping, because they are optimizing
different
>parameters. Substituting "z=exp(sT)" is a naive mapping of critical points,
>without regard to the effect upon magnitude response. Your method is an
explicit
>minimization of the differences between the prototype and model magnitude
>responses, without regard to the effects upon locations of poles and
zeroes.)
I am sure that there are other ways of doing this kind of thing, and for all
I know, one or more may be better in terms of stability, accuracy, or
computational efficiency. But my method has proven very useful to me because
it can produce results that are optimum in a minimax error sense over a
specified frequency range (and even if the results are not fully minimax
because of numerical issues, they are usually close), it is computationally
efficient, it can produce extremely accurate results with smooth curves like
RIAA and the like (often to within 1E-6 dB), and it almost always converges
when used to model classic analog magnitude responses required in audio.
For example, I have used it automatically inside a loop to create tables of
pole-zero values for parametric equalizers that model the magnitude responses
of classic analog second-order bell-curve equalizers. Each table contains
several hundred values that correspond to w0 and Q pairs and takes less than
a minute to compute with a reasonably modern PC. It is interesting to look at
the IIR coefficients in these tables because they are a bit chaotic. (A 3D
surface plot of a given IIR coefficient vs. frequency vs. Q isn't
particularly smooth.) I originally thought that it would be possible to do
some sort of linear interpolation between values, but it turns out that big
jumps occur. This is evidently a consequence of the minimax error criterion
but I didn't expect it before it saw it.
Here's an interesting example -- matching a 1 kHz 3rd order analog
Butterworth filter from 0 to 21 kHz at 44.1 kHz sample rate. (I used rounded
-off s-plane poles: -1, -0.5+/-j 0.866.) The target curve has an 80 dB
magnitude range, give or take. This is also challenging because the upper
frequency bound is 95.2% of Nyquist.
Here is the result of a 6-th order optimization:
Zero # Real Imag.
1 -0.8840140 0.000000
2 -0.7325760 0.000000
3 -0.5384182 0.000000
4 -0.3346849 0.000000
5 -0.5219714E-01 0.000000
6 0.3407240E-02 0.000000
Pole # Real Imag.
1 -0.8733325 0.000000
2 -0.6922994 0.000000
3 -0.4365723 0.000000
4 0.8672077 0.000000
5 0.9241607 0.1146088
6 0.9241607 -0.1146088
MAXIMUM ERROR FROM 0.00 Hz TO 21000.00 Hz IS 0.0000871dB
What is interesting here is what the algorithm does with the "extra" three z
-plane poles and zeros that it is given to work with. It seems to interlace
them on the real axis in the z-plane (see the first three poles and zeros
above). There is nothing explicit in the algorithm that forces this to
happen; it seems to happen naturally as part of the optimization process.
When the algorithm is allowed only three poles and zeros, it still gives good
results. It is interesting to compare the lower 3 poles and zeros in the 6-th
order optimization with the poles and zeros produced by the third-order
optimization:
Zero # Real Imag.
1 -0.6030164 0.000000
2 0.2847928E-01 0.2376939
3 0.2847928E-01 -0.2376939
Pole # Real Imag.
1 0.8486569 0.000000
2 0.9263526 0.1100756
3 0.9263526 -0.1100756
MAXIMUM ERROR FROM 0.00 Hz TO 21000.00 Hz IS 0.1393570dB
This one has complex zeros, which the 6-th order optimization did not.
A 4-th order optimization does not converge.
A 5-th order optimization converges but is not quite minimax because of ill-
conditioning in the Remez root finder. Nevertheless, it greatly improves the
error compared to the third-order solution:
Zero # Real Imag.
1 -0.8485930 0.000000
2 -0.6440151 0.000000
3 -0.3913983 0.000000
4 -0.1923411E-01 0.2816178E-01
5 -0.1923411E-01 -0.2816178E-01
Pole # Real Imag.
1 -0.8286350 0.000000
2 -0.5579574 0.000000
3 0.8671565 0.000000
4 0.9241636 0.1145946
5 0.9241636 -0.1145946
MAXIMUM ERROR FROM 0.00 Hz TO 21000.00 Hz IS 0.0004787dB
A 7-th order solution is the best the algorithm can find. It not quite
minimax because it encounters a root on the positive real axis (a forbidden
condition, as explained in my earlier post on the optimizer):
Zero # Real Imag.
1 -0.9072050 0.000000
2 -0.8028597 0.000000
3 -0.6672917 0.000000
4 -0.4965467 0.000000
5 -0.3171957 0.000000
6 -0.5621694E-01 0.000000
7 0.5741816E-02 0.000000
Pole # Real Imag.
1 -0.9011459 0.000000
2 -0.7853810 0.000000
3 -0.6211117 0.000000
4 -0.3974577 0.000000
5 0.8672086 0.000000
6 0.9241608 0.1146089
7 0.9241608 -0.1146089
MAXIMUM ERROR FROM 0.00 Hz TO 21000.00 Hz IS 0.0000147dB
The algorithm produces no improvement beyond 7-th order, but 0.0000147dB
error is good enough for most folks :-)
As far as common closed-form transformations go, the step-invariant transform
seems to be best:
Zero # Real Imag.
1 -3.473398 0.000000
2 -0.2496715 0.000000
Pole # Real Imag.
1 0.8672085 0.000000
2 0.9241609 0.1146089
3 0.9241609 -0.1146089
It provides +/- 1.1 dB accuracy over this frequency range. The result of the
direct calculation of the step-invariant transform is non-minimum phase but
could be easily corrected by reflecting the -3.473398 zero around the unit
circle. However, the zeros shown above provide the best phase match to the
analog prototype after a fixed delay is subtracted from the group delay to
match the group delay of the analog and digital filters as much as possible.
The phase match of the minimax solution to the analog filter is actually much
poorer compared to the step-invariant transform. Even though the minimax
approximation is minimum-phase (as is the analog filter), the phase response
is strongly affected by the magnitude response outside of the approximation
band. The group delay of the minimax digital filter becomes negative above
about 3.5 kHz. The group delay of the analog prototype is positive at all
frequencies, as is the group delay of the step-invariant transform. So this
is an interesting example of a case where a non-minimum phase digital
approximation matches the group delay of a minimum-phase analog filter better
than a minimum-phase approximation with the same magnitude response.
Matched-z and bilinear as basically useless, giving over 10 dB of error. The
next best is impulse-invariant (below), which gives about +/- 3.5 dB of
error.
Zero # Real Imag.
1 -0.9094103 0.000000
2 0.000000 0.000000
Pole # Real Imag.
1 0.8672085 0.000000
2 0.9241609 0.1146089
3 0.9241609 -0.1146089
--Bob Orban