[SciPy-dev] RE: Remez algorithm
Chuck Harris
Chuck.Harris at sdl.usu.edu
Wed Apr 24 10:49:35 CDT 2002
> Have you seen the remez algorithm in the signal module?
>
> How does this one compare?
>
> -Travis
>
This looks to be a pretty close implementation of the original
Parks-McClellan algorithm. As such it directly applies to real
symmetric filters. The idea is that in the frequency domain the
transfer function is a sum of terms of the form cos(2*pi*n*f)
which may be rearranged into powers of n, cos(2*pi*f)**n, so that
a change of variable to cos(2*pi*f) yields a Vandermonde matrix and
Lagrange interpolation can be used to evaluate sums of the form
a + b*cos(2*pi*1*f) + c*cos(2*pi*2*f) + ... passing thru a specified
set of points without ever having to solve for a,b,c, etc.
At the least, I feel that this should be generalized so as to accept
transfer functions like sum( a[n]*exp(2*pi*1j*n*f)), which are also
Vandermonde matrices after the change of variable x <- exp(2*pi*1j*f),
that is, it would be nice to use Lagrange interpolation for points
specified on the unit circle in the complex plane, which is why I sent
along nevilles version of Lagrange interpolation, which should perhaps
be redone in C for speed. So my first comment is that the routine should
be adapted to use a general, good form of Lagrange interpolation which
should be included in scipy.
The routine I sent *does* solve for the coefficients directly, and then
evaluates the sum by matrix multiplication. It is of slightly lower order
than using Lagrange interpolation when the number of points to be interpolated
at exceeds the number of coefficients, which is always, and is more general,
in that it is not necessary to have matrices that could be reduced to Vandermonde
form in principal. There is a generalization of Lagrange interpolation due to
Muhlbach that applies to complete Chebychev systems, and such might be an
interesting topic to pursue.
In any case, the linear algebra routines in SciPy are highly optimized and
are much faster than Lagrange interpolation in its present form. The
potential drawback is that the matrices appearing in the solution may
not be well conditioned, so some good sense is required on the part of
the user. One could argue that the badly conditioned case reflects a
certain indeterminacy in the coefficients that one wants to use, and
as such reflects reality. On the other hand, Lagrange interpolation will
produce pretty accurate interpolations whatever the users choice of basis,
and uses less space as no matrices have to be maintained.
The reason I wrote this routine is that Remez package available in MatLab
was also too restrictive and I needed something better.
Chuck
More information about the Scipy-dev
mailing list