Remez algorithm and filter lengths

Hi all. 

I was just involved in a thread over at comp.soft-sys.matlab where I got 
involved with the Remez algorithm. The matlab SP toolbox version is 
the only implementation of the Remez algorithm I have used, and it 
requires the user to specify the length of the FIR filter to be 
designed. 

Now, does the Remez algorithm generally require the filter length to 
be supplyed by the user, or does it exist versions that estimate 
the number of filter coefficients as well?

Rune

Hello Rune,
Specifically the Remez algorithm is one for finding the best fit polynomial
in a Chebyshev sense (I..e, minimize the maximum error). And yes the
polynomial order is prespecified. Now in the world of DSP, many call the
Parks-McClellan algo a Remez algo. But actually the PM algo is one that
takes a filter spec and represents as a poly in cosine functions, then uses
the Remez algo to find the best fit, and then converts the results back to
the domain of filters. I think of the PM algo as a wrapper that translates a
problem to one of poly approximation then uses Remez to find the poly and
then applies an inverse translation. I use the Remez algo for many things
other than filter problems. The basic idea behind PM is given below:

For example an N tap linear phase filter (even symmetry) will have have a
frequency repsonse in the functional form of a0 +
a1*cos(wt)+a2*cos(2wt)+a3*cos(3wt)+...

Now using some trig identities, we can rewrite this as
b0+b1*cos(wt)+b2*cos^2(wt)+b3*cos^3(wt)+...

Which after letting x=cos(wt), we get

b0+b1*x+b2*x^2+b3*x^3+...

which is simply a poly in x. So we use the Remeze algo to find the best fit
(which values of b0,b1,b2,...bn) make the poly come closest to the desired
freq. response in a min-max way.

Once these are found, the inverse of the original trig identities are used
to find the a0,a1,a2,... which simply relate to the filter taps. Of course
one can simply sample the best fit poly and do an inverse DFT. PM did this
in one of their early programs, but then they later used a computationally
more efficient method of inverting the trig identities.

For the other 3 cases of linear phase filters (they can be either symmetric
or antisymmetric and they can be either even or odd in length), a simple
manipulation can be done to put the approximation into one a fitting a
cosine poly.

The Hilbert transform and differentiators are sin functions and a single
trig term may be factored out leaving a cosine series. And the even length
cases are handled similarly.


But the answer to your question is one puts in the poly order (related to
the filter order - depends on even,odd and symmetry type) and then finds the
filter. But there are some equations (separate from PM) that allow you to
approximate the filter order given its specs. It is interesting to note that
sometimes an N-1 order filter may fit better than an Nth order one. This has
to do with the unification that puts all 4 types into one type of appx.
problem. And therefore the basis functions for the N-1 and N cases are
different.

Clay





"Rune Allnor" <allnor@tele.ntnu.no> wrote in message
news:f56893ae.0312221215.5f3f8f01@posting.google.com...
> Hi all.
>
> I was just involved in a thread over at comp.soft-sys.matlab where I got
> involved with the Remez algorithm. The matlab SP toolbox version is
> the only implementation of the Remez algorithm I have used, and it
> requires the user to specify the length of the FIR filter to be
> designed.
>
> Now, does the Remez algorithm generally require the filter length to
> be supplyed by the user, or does it exist versions that estimate
> the number of filter coefficients as well?
>
> Rune

"Rune Allnor" <allnor@tele.ntnu.no> wrote in message
news:f56893ae.0312221215.5f3f8f01@posting.google.com...
> Hi all.
>
> I was just involved in a thread over at comp.soft-sys.matlab where I got
> involved with the Remez algorithm. The matlab SP toolbox version is
> the only implementation of the Remez algorithm I have used, and it
> requires the user to specify the length of the FIR filter to be
> designed.
>
> Now, does the Remez algorithm generally require the filter length to
> be supplyed by the user, or does it exist versions that estimate
> the number of filter coefficients as well?

Rune,

Yes, the Remez exchange algorithm is a search algorithm that uses a fixed
number of terms.  What with the speed of the algorithm for most things, I've
often wondered why someone didn't package it up in an interative form that
would determine the length for a given set of specifications.  Seems it
would be easy enough:

1) Start with a Type as usual.  Use the estimate formulae for a starting
point.
2) if the result isn't sufficient/ is overkill, increase / decrease the
length to something longer / shorter that fits the Type.
Continue in a brute force way until the desired result is just achieved -
yielding the length as a by-product.

Hey, I'm tempted to modify the PM to do it.  But, before I do that, I'd like
to make sure whatever I end up with will be capable of generating very long
filters - because not all programs are.

Fred

In article f56893ae.0312221215.5f3f8f01@posting.google.com, Rune Allnor at
allnor@tele.ntnu.no wrote on 12/22/2003 15:15:

> Hi all. 
> 
> I was just involved in a thread over at comp.soft-sys.matlab where I got
> involved with the Remez algorithm. The matlab SP toolbox version is
> the only implementation of the Remez algorithm I have used, and it
> requires the user to specify the length of the FIR filter to be
> designed. 
> 
> Now, does the Remez algorithm generally require the filter length to
> be supplyed by the user, or does it exist versions that estimate
> the number of filter coefficients as well?

isn't there a MATLAB function called remezord()?

 REMEZORD  FIR order estimator (lowpass, highpass, bandpass, multiband)
    [N,Fo,Ao,W] = REMEZORD(F,A,DEV,Fs) finds the approximate order N,
    normalized frequency band edges Fo, frequency band magnitudes Ao and
    weights W to be used by the REMEZ function as follows:
        B = REMEZ(N,Fo,Ao,W)
    The resulting filter will approximately meet the specifications given
    by the input parameters F, A, and DEV.  F is a vector of cutoff
    frequencies in Hz, in ascending order between 0 and half the sampling
    frequency Fs. If you do not specify Fs, it defaults to 2.  A is a vector
    specifying the desired function's amplitude on the bands defined by F.
    The length of F is twice the length of A, minus 2 (it must therefore
    be even).  The first frequency band always starts at zero, and the last
    always ends at Fs/2.  It is not necessary to add these elements to the
    F vector.  DEV is a vector of maximum deviations or ripples allowable
    for each band.
 
    C = REMEZORD(F,A,DEV,FS,'cell') is a cell-array whose elements are the
    parameters to REMEZ.
 
    EXAMPLE
       Design a lowpass filter with a passband cutoff of 1500, a
       stopband cutoff of 2000Hz, passband ripple of 0.01, stopband ripple
       of 0.1, and a sampling frequency of 8000Hz:
 
          [n,fo,mo,w] = remezord( [1500 2000], [1 0], [0.01 0.1], 8000 );
          b = remez(n,fo,mo,w);
 
       This is equivalent to
          c = remezord( [1500 2000], [1 0], [0.01 0.1], 8000, 'cell');
          b = remez(c{:});
 
   CAUTION 1: The order N is often underestimated. If the filter does not
   meet the original specifications, a higher order such as N+1 or N+2 will.
   CAUTION 2: Results are inaccurate if cutoff frequencies are near zero
   frequency or the Nyquist frequency.
 
    See also REMEZ, KAISERORD.

r b-j

Excellent synopsis.  Thanks.

Robert

www.gldsp.com

"Clay S. Turner" <CSTurner@WSE.Biz> wrote:

>Hello Rune,
>Specifically the Remez algorithm is one for finding the best fit polynomial
>

( modify address for return email )

www.numbersusa.com
www.americanpatrol.com

<r_obert@REMOVE_THIS.hotmail.com> wrote in message
news:cd4huv4lcd6t0i1biojbn0dlgjlu0jncfq@4ax.com...
> Excellent synopsis.  Thanks.
>

Clay,

Succinct too!

It took me a while to figure out why PM bothered to go to the x^n form of
the polynomial because the cosine nx series is a perfectly good basis set.
It seems the barycentric form of Lagrange interpolation - which they use to
find the next approximant / design - needs to use that x^n form.  Otherwise,
they could have used the more direct form of cosine(nx) and generated the
filter coefficients directly.  I also don't understand why they used the -
what was it, 1/x? - weighting (I don't mean filter design criteria
weighting).

Anyway, I've always viewed the PM implementation with Remez exchange as the
"engine" as being pretty direct - but probably because I ignored the x^n
part inside - it was "just a detail".

In my understanding, the Remez exchange algorithm isn't limited to
polynomials.  Any sum of basis functions that satisfy the Haar condition
will do - and there are many.  Sets such as: cos(nx), sin(nx)/x, etc....

Also, there are many ways to solve for the next approximant:
1) Solution of n+1 linear equations in n coefficients and one error term.
2) Lagrange interpolation
3) Maybe some form of Newton interpolation - which can be more efficient in
some situations

I think the first one is hard to get to work on a PC with double precision
when the order gets to be in the 100's.
The second one may be better in this regard.
I don't know about the third one.

I'm interested in knowing what method people use to get long filters like
length 4096, etc.  Any suggestions?  I want to code up a better set of
programs.

Fred

On 22 Dec 2003 12:15:34 -0800, allnor@tele.ntnu.no (Rune Allnor)
wrote:

>Hi all. 
>
>I was just involved in a thread over at comp.soft-sys.matlab where I got 
>involved with the Remez algorithm. The matlab SP toolbox version is 
>the only implementation of the Remez algorithm I have used, and it 
>requires the user to specify the length of the FIR filter to be 
>designed. 
>
>Now, does the Remez algorithm generally require the filter length to 
>be supplyed by the user, or does it exist versions that estimate 
>the number of filter coefficients as well?
>
>Rune

Hi Rune,
   here's an example of what I use:

% Design lowpass FIR using remezord() & remez() commands. 
  Rp = 3;          % Desired passband ripple in dB
  Rs = 30;         % Desired stopband attenuation in dB
  Fs = 7000;       % Sampling frequency
  Freq = [400,500]; % End of passband & start of stopband
  Mag = [1,0];       % Desired magnitudes
% Compute deviations vector needed by the remezord() command.
  Dev = [(10^(Rp/20)-1)/(10^(Rp/20)+1),10^(-Rs/20)];
  [K,Freqo,Mago,Weight] = remezord(Freq,Mag,Dev,Fs)
  B = remez(K,Freqo,Mago,Weight);
  freqz(B,1,200), zoom on


[-Rick-]

Hello Fred,
Comments interspersed below:


"Fred Marshall" <fmarshallx@remove_the_x.acm.org> wrote in message
news:eomdnUxuDP4BCnWiRVn-hA@centurytel.net...
>
> <r_obert@REMOVE_THIS.hotmail.com> wrote in message
> news:cd4huv4lcd6t0i1biojbn0dlgjlu0jncfq@4ax.com...
> > Excellent synopsis.  Thanks.
> >
>
> Clay,
>
> Succinct too!

Sometimes I'm too lazy to write a lot, so my way is to be terse.

>
> It took me a while to figure out why PM bothered to go to the x^n form of
> the polynomial because the cosine nx series is a perfectly good basis set.
> It seems the barycentric form of Lagrange interpolation - which they use
to
> find the next approximant / design - needs to use that x^n form.
Otherwise,
> they could have used the more direct form of cosine(nx) and generated the
> filter coefficients directly.  I also don't understand why they used the -
> what was it, 1/x? - weighting (I don't mean filter design criteria
> weighting).

The weighting shows up with the type II,III, and IV filters. Basically you
have a sum of trig functions and part of the change to the polynomial from
involves the factoring out of a trig function. Since the poly sans the
factored out trig function is what is "best fitted", some compensation is
needed for the factored out trig function.

Examples:


Even Sym - Even Length

H(f) = cos(pi*F) * sum( b_k*cos(2pikF))   k is in 0,1,2,...,n-1

Odd Sym - Odd Length

H(f) = sin(2pi*F) * sum( b_k*cos(2pikF))   k is in 0,1,2,...,n-1

Odd Sym - Even Length

H(f) = sin(pi*F) * sum( b_k*cos(2pikF))   k is in 0,1,2,...,n-1

In all of the above cases, the cosine poly is what is "Remezed" - The
factored out term acts as a weight function.

>
> Anyway, I've always viewed the PM implementation with Remez exchange as
the
> "engine" as being pretty direct - but probably because I ignored the x^n
> part inside - it was "just a detail".
>
> In my understanding, the Remez exchange algorithm isn't limited to
> polynomials.  Any sum of basis functions that satisfy the Haar condition
> will do - and there are many.  Sets such as: cos(nx), sin(nx)/x, etc....

The Remez algo will find the best fit (in a min-max way) poly to fit the
data. That's how the algo is developed. However there is nothing to say you
can't change from one basis to another as long as certain criteria are
adhered to. Of course Remez himself proved that as you moved the abcissal
values about, the deviation factor rho takes on an extreme value when the
chebyshev criterion is met. His algo as he presented it does assume a
polynomial, but as PM showed, this polynomial can also represent something
else.

>
> Also, there are many ways to solve for the next approximant:
> 1) Solution of n+1 linear equations in n coefficients and one error term.
> 2) Lagrange interpolation
> 3) Maybe some form of Newton interpolation - which can be more efficient
in
> some situations
>
> I think the first one is hard to get to work on a PC with double precision
> when the order gets to be in the 100's.
> The second one may be better in this regard.
> I don't know about the third one.

1) Inverting a large matrix is quite painful and the numerical issues will
be tremendous

2) LaGrange Interpolation will work but the barycentric form is both more
efficient and better from a numerical point of view.

3) The Newton method is quite efficient but as you build a table adding in
successive points, the accuracy obtained in interpolation depends on whether
you are looking near the early used or the latter used points in the
tableau. This is not very good for high order interpolation. In Lagrange
interpolation, all of the points are equally represented in the
interpolation process, so the errors are spread out.

4) The Barycentric form is good in that it is a two step interpolation
method with the results from the 1st step being only dependent on the
abcissal values. Thus if you want to interpolate a bunch of points as is
done in the Remez algo, you do step 1 once, and then the step two operation
is done for each point. This actually results in less overall computation.
Plus there is a bunch of literature on the good numerical properties of
barycentric interpolation.

>
> I'm interested in knowing what method people use to get long filters like
> length 4096, etc.  Any suggestions?  I want to code up a better set of
> programs.

I've used double precision to design filters with length 2048. It would be a
simple matter to overload the math operators with a high precision math
library. I have my Remez written in c, but it would be zero effort to change
it over to c++ so it could use high precision math. With double precision
case, many failures stem from underflow[1]. I can't design filters with 200
dB attenuation. But if I keep the transition bandwiths narrow enough, I can
design very long filters.

I thought some designed the very big filters by windowing methods.

Clay

[1] The key calculations in interpolation involve differences and when these
throw away most of the digits (from using similarly sized numbers), this is
when the trouble arises.


>
> Fred
>
>
>
>
>
>