The Parks-McClellan algorithm

Hi All. 

To make a short story long, my little FIR filter project is happily 
trodding along. I am starting to get the various window functions in 
place for LP prototypes, and with some luck, it's not a big deal to 
generate HP or BP versions from the LP prototypes. Band-stop filters 
are a bit more cumbersome, I have to thing about how to implement them. 

All the FIR window design stuff appears to be going on well. So, I thougt
(infinitely naive as I am) why not have a go with the Parks-McClellan 
algorithm as well. 

I searched at IEEExplore and found four or five papers from the early 
1970ies authored by the venerable Parks & McClellan and found flow 
charts and FORTRAN code (spaghetti version) galore. I am not very fond of 
untangling spaghetti code, but the flow charts are clear enough for me to 
implement my own version of the algorithm. 

Except for the Remez exchange step that apparently is so crucial. So I 
consulted "Numerical Recipies" to see what they had to say about the 
Remez algorithm. They basically said "we don't like the Remez thing, 
we prefer to do things by brute force", and then they went on by 
explaining how to do things by brute force. 

To make the above story short, my question is if anyone knows of a 
readable and accesible account of the Remez exchange algorithm.

Rune

Rune Allnor wrote:
> To make the above story short, my question is if anyone knows of a 
> readable and accesible account of the Remez exchange algorithm.

The canonical comp.dsp reference for Remez would be Jake Janovetz's 
implementation:

http://www.janovetz.com/jake/

Jake used to be a regular here, and his work on Remez has received 
plenty of comp.dsp accolades.  I think his "account" of Remez is written 
in C.

-- 
Jim Thomas            Principal Applications Engineer  Bittware, Inc
jthomas@bittware.com  http://www.bittware.com    (603) 226-0404 x536
There's a fine line between clever and stupid

Rune Allnor wrote:

> Hi All. 
> 
> To make a short story long, my little FIR filter project is happily 
> trodding along. I am starting to get the various window functions in 
> place for LP prototypes, and with some luck, it's not a big deal to 
> generate HP or BP versions from the LP prototypes. Band-stop filters 
> are a bit more cumbersome, I have to thing about how to implement them. 
> 
> All the FIR window design stuff appears to be going on well. So, I thougt
> (infinitely naive as I am) why not have a go with the Parks-McClellan 
> algorithm as well. 
> 
> I searched at IEEExplore and found four or five papers from the early 
> 1970ies authored by the venerable Parks & McClellan and found flow 
> charts and FORTRAN code (spaghetti version) galore. I am not very fond of 
> untangling spaghetti code, but the flow charts are clear enough for me to 
> implement my own version of the algorithm. 
> 
> Except for the Remez exchange step that apparently is so crucial. So I 
> consulted "Numerical Recipies" to see what they had to say about the 
> Remez algorithm. They basically said "we don't like the Remez thing, 
> we prefer to do things by brute force", and then they went on by 
> explaining how to do things by brute force. 
> 
> To make the above story short, my question is if anyone knows of a 
> readable and accesible account of the Remez exchange algorithm.
> 
> Rune

"An Introduction to the Approximation of Functions" by Theodore J.
Rivlin. My Dover reprint cost $3.50, but the price has risen now to
$8.95. http://store.yahoo.com/doverpublications/0486640698.html

It's rather abstruse (Rick will look that up), maybe too hard for me,
but duck soup for you. What other book can be has for about the same
price as the postage?

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������

Rune Allnor wrote:
> Hi All. 
> 
> To make a short story long, my little FIR filter project is happily 
> trodding along. I am starting to get the various window functions in 
> place for LP prototypes, and with some luck, it's not a big deal to 
> generate HP or BP versions from the LP prototypes. Band-stop filters 
> are a bit more cumbersome, I have to thing about how to implement them. 
> 
> All the FIR window design stuff appears to be going on well. So, I thougt
> (infinitely naive as I am) why not have a go with the Parks-McClellan 
> algorithm as well. 
> 
> I searched at IEEExplore and found four or five papers from the early 
> 1970ies authored by the venerable Parks & McClellan and found flow 
> charts and FORTRAN code (spaghetti version) galore. I am not very fond of 
> untangling spaghetti code, but the flow charts are clear enough for me to 
> implement my own version of the algorithm. 
> 
> Except for the Remez exchange step that apparently is so crucial. So I 
> consulted "Numerical Recipies" to see what they had to say about the 
> Remez algorithm. They basically said "we don't like the Remez thing, 
> we prefer to do things by brute force", and then they went on by 
> explaining how to do things by brute force. 
> 
> To make the above story short, my question is if anyone knows of a 
> readable and accesible account of the Remez exchange algorithm.
> 
> Rune

An alternative is to use a general linear program (LP) solver.  There is 
a section in Rabiner and Gold that shows how, for Chebychev min max 
designs equivalent to Remez.  They, Rabiner and Gold, also make the 
statement that an advantage of the LP approach is being able to also 
include time domain constraints in the filter design.  The Remez 
approach is more efficient but since we don't think of of a VAX with 64K 
of memory as being a particularly representative workstation anymore, 
the differences in efficiency are not likely to be as critical as they 
once were.

Remez is a specialized algorithm. I don't know of any application of it 
  other than to Parks-Maclelan. Linear programming is considerably more 
general and INMHO more interesting and useful. Time spent understanding 
LP is well spent.  I'm not inclined to say the same of Remez.

There are also a number of LP solvers available on the web in various 
languages and in various shades of object orientation.

"Rune Allnor" <allnor@tele.ntnu.no> wrote in message
news:f56893ae.0409080215.48eff776@posting.google.com...
> Hi All.
>
> To make a short story long, my little FIR filter project is happily
> trodding along. I am starting to get the various window functions in
> place for LP prototypes, and with some luck, it's not a big deal to
> generate HP or BP versions from the LP prototypes. Band-stop filters
> are a bit more cumbersome, I have to thing about how to implement them.
>
> All the FIR window design stuff appears to be going on well. So, I thougt
> (infinitely naive as I am) why not have a go with the Parks-McClellan
> algorithm as well.
>
> I searched at IEEExplore and found four or five papers from the early
> 1970ies authored by the venerable Parks & McClellan and found flow
> charts and FORTRAN code (spaghetti version) galore. I am not very fond of
> untangling spaghetti code, but the flow charts are clear enough for me to
> implement my own version of the algorithm.
>
> Except for the Remez exchange step that apparently is so crucial. So I
> consulted "Numerical Recipies" to see what they had to say about the
> Remez algorithm. They basically said "we don't like the Remez thing,
> we prefer to do things by brute force", and then they went on by
> explaining how to do things by brute force.
>
> To make the above story short, my question is if anyone knows of a
> readable and accesible account of the Remez exchange algorithm.
>
> Rune

Hello Rune,

Here is "c" version of the Remez algo from the PM papers. Most of the
variable names are the same so deciphering it shouldn't be too hard.

The function Ak() calculates the Barycentric weights.

The functions BF() and BFc() evaluate the current LaGrange polynomial. BF()
is past an index to the grid which refers to the x value. Whereas BFc() is
simply passed an x value.

rhofun() calculates the chebyshev deviation factor.

I'm sorry the documentation is nonexistant. I haven't done anything with
this code in nearly 20 years. I did put together a complete PM algo. And I
can make that code available if desired.

Clay




/* ============================================================ */
void makeAk(int n)
{
int i,k;
long double q;
for (k=0;k<n;k++)
 if (k&1)
  ak[k]= -1.0;
 else
  ak[k]=  1.0;
for (k=0;k<n-1;k++)
 for (i=k+1;i<n;i++) {
  q=X[k]-X[i];
  ak[k]/=q;
  ak[i]/=q;
  }
}
/* ============================================================ */
DATA BF(int ind,int lg) /* barycentric lagrangian interpolation */
{
long double a,t,b,x1;
int k;
for (k=0;k<lg;k++)
 if (iext2[k]==ind)
  return(Y[k]);
x1=gridx[ind];
for (k=0,t=0.0,b=0.0;k<lg;k++) {
 a= ak[k]/(x1-X[k]);
 t+=(Y[k]*a);
 b+=a;
 }
return(t/b);
}
/* ============================================================ */
DATA BFc(DATA x,int lg)
{
long double a,t,b;
int k;
for (k=0;k<lg;k++)
 if (fabs(x-X[k])<0.000000001)
  return(Y[k]);
for (k=0,t=0.0,b=0.0;k<lg;k++) {
 a=ak[k]/(x-X[k]);
 t+=(Y[k]*a);
 b+=a;
 }
return(t/b);
}
/* ============================================================ */
DATA rhofun()
{
int j,k;
DATA top,bot;
 /* make Ak for rho function */
 makeAk(order+2);
 /* make the rho function */
 for (k=0,top=0.0,bot=0.0;k<order+2;k++) {
  j=iext[k];
  top+=(ak[k]*gridy[j]);
  if (k&1)
   bot-=(ak[k]/gridw[j]);
  else
   bot+=(ak[k]/gridw[j]);
  }
/*printf("top= %e   bot= %e    rho=  %1.16f order =
%d\n",top,bot,top/bot,order);*/
return(top/bot);
}
/* ============================================================ */
void remez(void)
{
long double x,y,err,f1,f2,meps=0.00000001;
int j,k,jlow,jup,jmax,iter=0,jchng=1,expts,numpts;
int forw,back,lg,flg1,flg2,ok;
numpts=order+2;
lg=order+1;
jmax=ngrid-1;


/* make order+2 initial guesses at the extremal points */
x=((long double)(jmax-2)/(long double)(lg)); /* numpts should = order+2  */
for (j=0,y=0.0;j<numpts;j++) {
 iext[j]=(int)(0.5+y);
 y+=x;
 }

/* do the remez approximation */

while (jchng && iter<17) {
 jchng=0;

 for (j=0;j<numpts;j++)
  iext2[j]=iext[j];

  /* snake polynomial through extremal points */
 for (j=0;j<numpts;j++)
  X[j]=gridx[iext[j]];
 rho=rhofun();
 for (j=0;j<numpts;j++) {
  k=iext[j];
  if (j&1)
   Y[j]=gridy[k]+rho/gridw[k];
  else
   Y[j]=gridy[k]-rho/gridw[k];
  }

 makeAk(lg); /* make LaGrangian interpolation coefficients */

 for (expts=0;expts<numpts;expts++) {


  /* set up search limits for current point */
  if (!expts) {
   jlow=0;
   jup=iext[1]-1;
   }
  else if (expts==lg) {
   jlow=iext[order]+1;
   jup=jmax;
   }
  else {
   jlow=iext[expts-1]+1;
   jup=iext[expts+1]-1;
   }
  j=iext[expts]; /* set up current pointer */
  forw=0;
  back=0;
  if (j==jlow)
   forw++; /* indicate a forced forward search */
  if (j==jup)
   back++; /* indicate a forced backward search */
  err= (Y[expts]-gridy[j])*gridw[j];

  if ((err>0.0) && (jup!=jlow)) {  /* looking for maximum */
   if (! (forw || back)) /* is there a forced direction? */
    if ( (BF(j+1,lg)-gridy[j+1])*gridw[j+1]-err>meps)
     forw++;
   if (forw) {  /* search forwards */
    j++; /* advance to interpolated point */
    ok=1;
    while (j<=jup && ok)
     if ((y=(BF(j,lg)-gridy[j])*gridw[j])-err>meps) {
      j++;
      err=y;
      jchng++;
      }
     else
      ok=0;
    j--;
    }

   else {  /* search backwards */
    j--; /* backup to interpolated point */
    ok=1;
    while (j>=jlow && ok)
     if ((y=(BF(j,lg)-gridy[j])*gridw[j])-err>meps) {
      j--;
      err=y;
      jchng++;
      }
     else
      ok=0;
    j++;
    }
   }
  else if (jup!=jlow) {  /* looking for minimum */
   if (! (forw || back)) /* is there a forced direction? */
    if (err-(BF(j+1,lg)-gridy[j+1])*gridw[j+1]>meps)
     forw++;
   if (forw) {  /* search forwards */
    ok=1;
    j++; /* advance to interpolated point */
    while (j<=jup && ok)
     if (err-(y=(BF(j,lg)-gridy[j])*gridw[j])>meps) {
      j++;
      err=y;
      jchng++;
      }
     else
      ok=0;
    j--;
    }

   else {  /* search backwards */
    j--; /* backup to interpolated point */
    ok=1;
    while (j>=jlow && ok)
     if (err-(y=(BF(j,lg)-gridy[j])*gridw[j])>meps) {
      j--;
      err=y;
      jchng++;
      }
     else
      ok=0;
    j++;
    }
   }
  iext[expts]=j; /* save the newly found point */

  } /* end of point search */
 flg1=0;
 if(iext[0]==0)
  flg1=1;
 flg2=0;
 if (iext[lg]==jmax)
  flg2=1;

 if (flg1^flg2) { /* enter iff one end point is max */
  if (flg2) {  /* right end is max */
   x=gridw[jmax]*(BF(jmax,lg)-gridy[jmax]);
   y=gridw[0]*(BF(0,lg)-gridy[0]);
   f1=BF(iext[0],lg)-gridy[iext[0]];
   if ( fabs(y)>1.0000001*fabs(x) ) {
    if (y*f1<0.0) {
     for (j=lg;j>0;j--)
      iext[j]=iext[j-1];
     iext[0]=0;
     jchng++;
     }
    }
   }
  else {   /* left end is max */
   x=gridw[0]*(BF(0,lg)-gridy[0]);
   y=gridw[jmax]*(BF(jmax,lg)-gridy[jmax]);
   f2=BF(iext[lg],lg)-gridy[iext[lg]];
   if (fabs(y)>1.0000001*fabs(x))
    if (y*f2<0.0) {
     for (j=0;j<lg;j++)
      iext[j]=iext[j+1];
     iext[lg]=jmax;
     jchng++;
     }
   }
  }

 else if (!(flg1|flg2)) { /* come here if neither end point is maximum */
  x=gridw[0]*(BF(0,lg)-gridy[0]);
  f1=BF(iext[0],lg)-gridy[iext[0]];
  y=gridw[jmax]*(BF(jmax,lg)-gridy[jmax]);
  f2=BF(iext[lg],lg)-gridy[iext[lg]];
  if (fabs(y)>fabs(x)) {
   if (y*f2<0.0) {
    for (j=0;j<lg;j++)
     iext[j]=iext[j+1];
    iext[j]=jmax;
    jchng++;
    }
   }
  else
   if (x*f1<0.0) {
    for (j=lg;j>0;j--)
     iext[j]=iext[j-1];
    iext[j]=0;
    jchng++;
    }
  }
 iter++;
 } /* end of remez approximation */
}

"Rune Allnor" <allnor@tele.ntnu.no> wrote in message
news:f56893ae.0409080215.48eff776@posting.google.com...
> Hi All.
>
> To make a short story long, my little FIR filter project is happily
> trodding along. I am starting to get the various window functions in
> place for LP prototypes, and with some luck, it's not a big deal to
> generate HP or BP versions from the LP prototypes. Band-stop filters
> are a bit more cumbersome, I have to thing about how to implement them.
>
> All the FIR window design stuff appears to be going on well. So, I thougt
> (infinitely naive as I am) why not have a go with the Parks-McClellan
> algorithm as well.
>
> I searched at IEEExplore and found four or five papers from the early
> 1970ies authored by the venerable Parks & McClellan and found flow
> charts and FORTRAN code (spaghetti version) galore. I am not very fond of
> untangling spaghetti code, but the flow charts are clear enough for me to
> implement my own version of the algorithm.
>
> Except for the Remez exchange step that apparently is so crucial. So I
> consulted "Numerical Recipies" to see what they had to say about the
> Remez algorithm. They basically said "we don't like the Remez thing,
> we prefer to do things by brute force", and then they went on by
> explaining how to do things by brute force.
>
> To make the above story short, my question is if anyone knows of a
> readable and accesible account of the Remez exchange algorithm.

Rune,

You might take a look at that IEEE paper reference I just gave you.  I think
the step by step description of the algorithm was pretty easy to read!
Then, if you like, you can replace the linear equation solution step with
some other method.

Here's what it says:  N is the number of filter coefficients

1) Start with a set of N+1 points - often equidistant over the range of
approximation.

2) Solve a set of linear equations in the filter coefficients *and* the
equalized magnitude of the error as the variables with the sign of the error
set up to alternate as required.  This results in an approximant with error
of sequentially alternating signs.  (The "approximant" is the current filter
frequency response).  You can introduce weighting of the error in these
equations if you wish - high weight, low absolute error - usually dependent
on the location of the N+1 points.

3) Using this approximant, compute the new error on a fine grid - relative
to the ideal filter response that you want ...  just a simple subtraction.

4) Find a new set of error maxima on the fine grid that meet the sign
alternation criterion (see below).

5) Using the indices of the peaks obtained in step 4, repeat the process
from step 2 until the error peaks in step 4 are all equal (within some small
range).  This is the unique minimax solution - which may have more than N+1
equal peaks.

I will say this:
Parks & McClellan used the barycentric form of Lagrange interpolation to
arrive at the next approximant.  I used a standard linear equation solver.

A lot of the fun is how you select the next points for the following
computation which will equalize the error on those points.

The perhaps simplest and slowest method (at least I imagine it must be the
slowest) is to simply exchange in the new highest peak for the adjacent peak
of the same sign - keeping N of the original N+1 points.

A faster method I believe (in terms of outer loop iteration count) is to
find *the* set of highest peaks that meets the sign alternation criterion.

If you're going to build your own Remez exchange algorithm program - that's
where the work will be.  Finding the peaks of the error is pretty easy.  I
did it in sort of a brute force way that I found satisfying - but perhaps
without tremendous computational efficiency.  (Parks & McClellan were more
clever and took advantage of some known attributes of the functions being
dealt with - mine is more brute force).

1) find *all* the points where the slope of the error changes sign (the
"peaks")

2) remove all of those from the sequence above that are contiguous and have
the same sign of the error - keeping only the one in each such sequence with
the highest magnitude.
This guarantees that you have at least the number of peaks necessary - and
that you have the largest ones.

3) As long as the number of remaining peaks is N+3 or greater, remove the
smallest peak from the remaining sequence.  Doing this will then force you
to remove the smallest adjacent peak of the opposite sign - to keep the sign
alternation property. You cannot remove the largest peak of the entire
sequence by doing this - which is important.

4) If you get to N+2 peaks, remove the smallest peak from one end of the
sequence.

This gives you *the* set of N+1 largest peaks that meet the sign alternation
criterion.  (Well, I say *the* set ... it may not be unique e.g. it could
happen that extra end points have the same magnitude).

I hope this helps.

Fred

And, I might have mentioned something that Parks and McClellan knew and it
took me a while to figure out:

You don't have to include the "don't care" regions (e.g transition bands) in
the region of approximation.  You can consider the discontinuous but
contiguous stop and pass bands as a continuum for purposes of finding the
error peaks.  Somehow this works out.  Maybe it's in Rabiner's book ....  At
least when I've tried it, it seems to work and reading the P-M code that's
how they implement it I do believe.

Alternately one could include the transition bands with a very low weight
but I find that this can lead to peaks occurring in the transition bands -
which is inconvenient because they don't belong there.  It all has to do
with the alternation of signs in the error.

It's easy enough to envision:  the passband amplitude goes into the
transition region with a value that's below ideal and thereafter enters the
stopband with a value that's above ideal.  If all of the peaks are in the
pass and stop bands, then there can't be any error peaks in the transition
bands...... something like that.

Fred

Hello Rune,

I put a .zip file of an application with the Remez algo that finds the best
poly fit to a user defined function. The code is a visual c++ project.
Go to the following link and download PolyApprox.zip.

http://personal.atl.bellsouth.net/p/h/physics/


Have fun,
Clay S. Turner

allnor@tele.ntnu.no (Rune Allnor) wrote in message news:<f56893ae.0409080215.48eff776@posting.google.com>...
> Hi All. 
> 
> To make a short story long, my little FIR filter project is happily 
> trodding along. I am starting to get the various window functions in 
> place for LP prototypes, and with some luck, it's not a big deal to 
> generate HP or BP versions from the LP prototypes. Band-stop filters 
> are a bit more cumbersome, I have to thing about how to implement them. 
> 
> All the FIR window design stuff appears to be going on well. So, I thougt
> (infinitely naive as I am) why not have a go with the Parks-McClellan 
> algorithm as well. 
> 
> I searched at IEEExplore and found four or five papers from the early 
> 1970ies authored by the venerable Parks & McClellan and found flow 
> charts and FORTRAN code (spaghetti version) galore. I am not very fond of 
> untangling spaghetti code, but the flow charts are clear enough for me to 
> implement my own version of the algorithm. 
> 
> Except for the Remez exchange step that apparently is so crucial. So I 
> consulted "Numerical Recipies" to see what they had to say about the 
> Remez algorithm. They basically said "we don't like the Remez thing, 
> we prefer to do things by brute force", and then they went on by 
> explaining how to do things by brute force. 
> 
> To make the above story short, my question is if anyone knows of a 
> readable and accesible account of the Remez exchange algorithm.
> 
> Rune

Try:
http://www.janovetz.com/jake/

SP [Wed, 08 Sep 2004 10:53:36 -0400]:
 >There are also a number of LP solvers available on the web in various 
 >languages and in various shades of object orientation.

Here's my favorite.  Man is that old, but it's pretty simple to
operate.  Does mixed and all-integer, too.  Click on the pic for
an example output.

 http://40th.com/dos/lp26.html

It comes with a nice little lpdemo.exe to get anyone started doing LP.
-- 
 40th Floor - Software  @  http://40th.com/
 iPlay : the ultimate audio player for iPAQs
 mp3, ogg, mp4, m4a, aac, wav, and then some