Remez algorithm and filter lengths

Started by Rune Allnor December 22, 2003
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 > > > > > >