Reply by Bogdan Kosanovic March 5, 20042004-03-05
You should try using Optimization toolbox in MATLAB. You need to define your
filter mask as a set of constraints (upper/lower masks) parametrized with
frequencies of interest. Your FIR coefficients should be optimized to
satisfy the desired response (e.g. middle of upper/lower mask) in Least
Square sense given the constraints that

|H(e^(jw))| <= upper mask-epsilon
-|H(h^(jw)| <= -(lower mask+epsilon)
|h(n)| < 1

Your "cost" function is J=|H(e^(jw)) - middle_line|^2

You may also add a gain parameter in case you need to move things up/down.
If you plan to use floating point implementation, no need for gain and
probably no need to constrain the h(n)'s to any specific range.

I guess fmincon() is the function to use from Optimization toolbox.

You may look deeper in Optimization toolbox and you could find a way to
apply arbitrary weights to each frequency of interest for each constraint.

If you want to add linear phase constraint, all you need to do is allow only
half of the coefficients to be optimized and the other half would be a
symmetric copy of those they you optimize.

Regards,
Bogdan

"Paul Mennen" <nospam@mennen.org> wrote in message
news:inw0c.18609$TJ5.4369@newssvr29.news.prodigy.com...
> Trying to design an FIR Low pass filter. > As a simple example, consider: > > >> edge = [0 .3 .7 1]; > >> weight = [100 1]; > >> 1000*remez(7,edge,[1 1 0 0],weight) > ans= 32.2690 -109.2208 81.9170 495.0689 495.0689 81.9170 -109.2208 32.2690 > > The resulting filter has a pass band ripple of +/- .013dB > and a stop band of -16.84dB. (The pass band is from DC to .3 > and the stop band is from .7 to 1 where 1 represents Fs/2). > > Note that the impulse response is symmetric so this is a linear phase > filter. This is expected because that is what REMEZ advertises itself > as doing. > > However I know I can get better performance by removing the linear > phase restriction. So now I try: > > >> 1000*cremez(7,edge,'lowpass',weight,'real') > ans= 32.2466 -109.1672 81.8522 495.0968 495.0968 81.8522 -109.1672 32.2466 > > >> 1000*cremez(7,edge,'lowpass',weight,'even') > ans= 32.2466 -109.1672 81.8522 495.0968 495.0968 81.8522 -109.1672 32.2466 > > Note that both of these responses are just about exactly the > same answer as given by REMEZ. > > From the help message (help cremez): > REMEZ Complex and nonlinear phase equiripple FIR filter design. > > Yet I see only linear phase filters coming out of cremez. > Isn't this false advertising? (I tried 'none' for the symmetry > argument, but a complex filter comes out - clearly not what > I am looking for.) > > The reason that I know that one can do much better than cremez > is that by trial and error I came up with a filter with this > impulse response: > > p = [307 535 307 -92 -113 43 29 -16]; > > I plotted the result with: plot(20*log10(abs(freqz(p,1,1000)))); > > Note the passband ripple is +/-.0126 dB (slightly better than before) > and the stopband is -24 dB (much better than before!). Clearly > not linear phase since it is not symmetrical, but I don't care > about that. > > Any ideas why my trial and error approach does so much better > than cremez? > > Any ideas on how to design such FIR filters without trial and error? > I thought of using yulewalk, but that has some problems. Yulewalk > is for iir filters, so right there it is at least twice as complicated > as it should be. Also Yulewalk minimizes in the least squares sense > and so it won't necessarily minimize the max deviation from desired. > And finally Yulewalk doesn't allow a weighting function. > > ~Paul Mennen >
Reply by Paul Mennen March 3, 20042004-03-03
> "Ricardo Losada" <Ric.Losada@mathworks.com> wrote > I want to address your post in two parts. ...
> Let me summarize first and then go into a little more detail.
Thanks Ricardo for the knowledgeable and informative reply! I agree that gremez looks like the way to go. I'll check it out. Too bad gremez is not included in the "see also" function list at the end of remez and cremez help. I suppose that was not done because gremez is in an optional toolbox, although it would certainly have saved me some time if I had been aware of that function sooner. Thanks again. ~Paul Mennen
Reply by Ricardo Losada March 2, 20042004-03-02
Hi,

I want to address your post in two parts. Let me summarize first and 
then go into a little more detail.

Bottom line:
CREMEZ can be used to design nonlinear phase filters as stated in the 
description of the function. However, it will probably not be useful for 
what you want.

On the other hand, your expectation of getting better performance from 
an FIR filter if you remove the linear-phase requirement is correct and 
there are tools from The MathWorks -specifically in the Filter Design 
Toolbox- to address this.

Details:
First let me address the concern on CREMEZ not being able to design 
nonlinear phase filters.

This can indeed be done, what needs to be done is use the delay 
parameter in the response to reduce the group delay compared to a linear 
phase design.

In your example:
b1 = 1000*cremez(7,edge,'lowpass',weight);

Which produces, as you mentioned, a linear phase filter with group delay 
of 3.5 samples.

If you specify a negative delay parameter (an advance effectively) like:
b2 = 1000*cremez(7,edge,{'lowpass',-2},weight);

The group delay is clearly not flat, thus the phase is nonlinear. The 
group delay tends to be about 1.5 (3.5+(-2)) in the passband of the 
filter. The impulse response is obviously no longer symmetric.

Unfortunately, this reduced passband group delay comes at the expense of 
decreased passband/stopband performance in the magnitude response. Which 
is the opposite of what you are looking for if I understand correctly.

What you are looking for can be achieved by using design functions 
available in the Filter Design Toolbox.

For instance, if you allow for even-ordered filters, if we compare these 
2 designs:
b1 = remez(8,edge,[1 1 0 0],weight);
b2 = gremez(8,edge,[1 1 0 0],weight,'minphase');

The nonlinear phase (minimum-phase) filter b2 has a better stopband 
attenuation (about 22 dB as opposed to just over 16 db) and less 
passband peak ripple (about 0.0014 dB compared to about 0.013 dB).

This is probably the way to go for equiripple designs if eve-ordered 
filters are ok.

If you absolutely need odd-order, or if you want to optimize for another 
norm, like a least-squares optimization, or a 4-norm, 8-norm, etc, it is 
best to use FIRLPNORM. For example:
b1=remez(7,edge,[1 1 0 0],weight);
b2=firlpnorm(7,edge,edge,[1 1 0 0],[100 100 1 1]);

In this case b2 has a stopband attenuation of about 24.48 dB as opposed 
to 16.87dB and a passband peak ripple of about 0.0052 dB as opposed to 
0.013 dB.

Finally, for a least squares response. If you compare:
b1=firls(7,edge,[1 1 0 0],weight);
b2=firlpnorm(7,edge,edge,[1 1 0 0],[10 10 1 1],[2 2]);

b2 will definitively have a smaller mean squared error than b1.

Hope this helps,
Ricardo.


Paul Mennen wrote:
> Trying to design an FIR Low pass filter. > As a simple example, consider: > > >>>edge = [0 .3 .7 1]; >>>weight = [100 1]; >>>1000*remez(7,edge,[1 1 0 0],weight) > > ans= 32.2690 -109.2208 81.9170 495.0689 495.0689 81.9170 -109.2208 32.2690 > > The resulting filter has a pass band ripple of +/- .013dB > and a stop band of -16.84dB. (The pass band is from DC to .3 > and the stop band is from .7 to 1 where 1 represents Fs/2). > > Note that the impulse response is symmetric so this is a linear phase > filter. This is expected because that is what REMEZ advertises itself > as doing. > > However I know I can get better performance by removing the linear > phase restriction. So now I try: > > >>>1000*cremez(7,edge,'lowpass',weight,'real') > > ans= 32.2466 -109.1672 81.8522 495.0968 495.0968 81.8522 -109.1672 32.2466 > > >>>1000*cremez(7,edge,'lowpass',weight,'even') > > ans= 32.2466 -109.1672 81.8522 495.0968 495.0968 81.8522 -109.1672 32.2466 > > Note that both of these responses are just about exactly the > same answer as given by REMEZ. > > From the help message (help cremez): > REMEZ Complex and nonlinear phase equiripple FIR filter design. > > Yet I see only linear phase filters coming out of cremez. > Isn't this false advertising? (I tried 'none' for the symmetry > argument, but a complex filter comes out - clearly not what > I am looking for.) > > The reason that I know that one can do much better than cremez > is that by trial and error I came up with a filter with this > impulse response: > > p = [307 535 307 -92 -113 43 29 -16]; > > I plotted the result with: plot(20*log10(abs(freqz(p,1,1000)))); > > Note the passband ripple is +/-.0126 dB (slightly better than before) > and the stopband is -24 dB (much better than before!). Clearly > not linear phase since it is not symmetrical, but I don't care > about that. > > Any ideas why my trial and error approach does so much better > than cremez? > > Any ideas on how to design such FIR filters without trial and error? > I thought of using yulewalk, but that has some problems. Yulewalk > is for iir filters, so right there it is at least twice as complicated > as it should be. Also Yulewalk minimizes in the least squares sense > and so it won't necessarily minimize the max deviation from desired. > And finally Yulewalk doesn't allow a weighting function. > > ~Paul Mennen >
Reply by Paul Mennen March 2, 20042004-03-02
"Fred Marshall" wrote: 
> First of all, you have to figure out if the linear phase > character of the filter is *really* costing you very much. > Maybe it isn't.
Gee I thought the simple example I gave in my original post showed pretty clearly that the cost of restricting a filter to linear phase is substantial. I've found that in higher order cases when one starts running into the limitations of the precision of the coefficients due to the number representation, this "cost" is even more dramatic.
> 1) Start by designing a linear phase FIR filter that ... > The method of defining an all-positive filter response and > taking the square root is due to Schuessler.
Thanks for describing this procedure. I had seen something like this years ago, (probably the same thing) although I had forgotten too many of the details to make use of it. I'm not sure how well the Schuessler procedure will work for my current problem since the design of the filter must be automated and this procedure might prove difficult to automate. That's why I was hoping the answer would just drop out of something like cremez (which seemingly is designed for just this sort of thing). ~Paul
Reply by Fred Marshall March 2, 20042004-03-02
"robert bristow-johnson" <rbj@surfglobal.net> wrote in message
news:BC698A2A.9104%rbj@surfglobal.net...
> In article NLqdnYxyQdCAcN7dRVn-uw@centurytel.net, Fred Marshall at > fmarshallx@remove_the_x.acm.org wrote on 03/01/2004 21:08: > > > > OK - so it's a matter of RTFM? > > dunno what the acronym is for.
R B-J, ***A software community term I believe: RTFM => "Read the f-ing manual" :-) So, I meant "It's just a matter of reading about "remez" and "firls" isn't it?" (referring to the OP). Sorry if that seems unkind. I am trying to be helpful. Fred
Reply by robert bristow-johnson March 2, 20042004-03-02
In article NLqdnYxyQdCAcN7dRVn-uw@centurytel.net, Fred Marshall at
fmarshallx@remove_the_x.acm.org wrote on 03/01/2004 21:08:

> > "robert bristow-johnson" <rbj@surfglobal.net> wrote in message > news:BC691892.90BD%rbj@surfglobal.net... >>> "A. Ser." wrote in message news:4043376A.8CD98C21@i.am... >>> >>> ..................... >>>> And here the algorithms like REMEZ >>>> often can give very bad results >>>> because it optimizes not the max. error >>>> but the mean square error. >>> ...................... >> >> that is not an accurate statement. >> >> In article Q4udncIskZWBPN7dRVn-tw@centurytel.net, Fred Marshall at >> fmarshallx@remove_the_x.acm.org wrote on 03/01/2004 15:46: >> >>> Well, I suppose that one could label *any* algorithm improperly. > However, >>> the Remez exchange algorithm *is* directed at minimizing the maximum > error. >>> Now, I can't vouch for what may be in Matlab or any other "collection" > with >>> such labels attached. >> >> well, i can. the MATLAB function "remez()" (Sig. Proc. Toolbox) > implements >> what would be more accurately called the "Parks-McClellan Algorithm" which >> uses the Remes Exchange Algorithm (which exists in Approximation Theory >> outside of any concept of filtering or DSP) to minimize the maximum > weighted >> error in designing FIR coefficients. the MATLAB function "firls()" does > the >> same thing (it is invoked exactly the same way) except that it minimizes > the >> total (or mean) square error. > > R B-J, > > OK - so it's a matter of RTFM?
dunno what the acronym is for.
> Remez for PM is unfortunate but probably not > a killer since that's what's buried inside. In my view PM did a very nice > job of wrapping the Remez exchange algorithm with some very nice practical > filter design interfaces and implementations - so it depends a lot on how > the interface is presented. FIRLS seems descriptive enough. So, my comment > to the OP remains....
fine, i didn't take issue with what you said. i *did* take issue with what A. Ser. said in that MATLAB's remez minimizes mean-square error. dunno who the OP is. r b-j
Reply by Fred Marshall March 1, 20042004-03-01
"robert bristow-johnson" <rbj@surfglobal.net> wrote in message
news:BC691892.90BD%rbj@surfglobal.net...
> > "A. Ser." wrote in message news:4043376A.8CD98C21@i.am... > > > > ..................... > >> And here the algorithms like REMEZ > >> often can give very bad results > >> because it optimizes not the max. error > >> but the mean square error. > > ...................... > > that is not an accurate statement. > > In article Q4udncIskZWBPN7dRVn-tw@centurytel.net, Fred Marshall at > fmarshallx@remove_the_x.acm.org wrote on 03/01/2004 15:46: > > > Well, I suppose that one could label *any* algorithm improperly.
However,
> > the Remez exchange algorithm *is* directed at minimizing the maximum
error.
> > Now, I can't vouch for what may be in Matlab or any other "collection"
with
> > such labels attached. > > well, i can. the MATLAB function "remez()" (Sig. Proc. Toolbox)
implements
> what would be more accurately called the "Parks-McClellan Algorithm" which > uses the Remes Exchange Algorithm (which exists in Approximation Theory > outside of any concept of filtering or DSP) to minimize the maximum
weighted
> error in designing FIR coefficients. the MATLAB function "firls()" does
the
> same thing (it is invoked exactly the same way) except that it minimizes
the
> total (or mean) square error.
R B-J, OK - so it's a matter of RTFM? Remez for PM is unfortunate but probably not a killer since that's what's buried inside. In my view PM did a very nice job of wrapping the Remez exchange algorithm with some very nice practical filter design interfaces and implementations - so it depends a lot on how the interface is presented. FIRLS seems descriptive enough. So, my comment to the OP remains.... Fred
Reply by robert bristow-johnson March 1, 20042004-03-01
> "A. Ser." wrote in message news:4043376A.8CD98C21@i.am... > > ..................... >> And here the algorithms like REMEZ >> often can give very bad results >> because it optimizes not the max. error >> but the mean square error. > ......................
that is not an accurate statement. In article Q4udncIskZWBPN7dRVn-tw@centurytel.net, Fred Marshall at fmarshallx@remove_the_x.acm.org wrote on 03/01/2004 15:46:
> Well, I suppose that one could label *any* algorithm improperly. However, > the Remez exchange algorithm *is* directed at minimizing the maximum error. > Now, I can't vouch for what may be in Matlab or any other "collection" with > such labels attached.
well, i can. the MATLAB function "remez()" (Sig. Proc. Toolbox) implements what would be more accurately called the "Parks-McClellan Algorithm" which uses the Remes Exchange Algorithm (which exists in Approximation Theory outside of any concept of filtering or DSP) to minimize the maximum weighted error in designing FIR coefficients. the MATLAB function "firls()" does the same thing (it is invoked exactly the same way) except that it minimizes the total (or mean) square error. people here know that i do not sing the praises of MATLAB because of their asinine hard-wired "1-based" indexing that allows no one, including God, to change that. but they are not false advertising that remez() uses the minimax norm to define error. r b-j
Reply by Fred Marshall March 1, 20042004-03-01
"Paul Mennen" <nospam@mennen.org> wrote in message
news:Csy0c.18642$lk7.4320@newssvr29.news.prodigy.com...

.............. That is why I was trying
> to figure out how to design non-linear phase FIR filters.
First of all, you have to figure out if the linear phase character of the filter is *really* costing you very much. Maybe it isn't. I guess the way to figure it out is to design a FIR filter that isn't linear phase. A good selection would be a minimum phase filter for lack of any other criterion. Now, there's a trivial way to go from linear phase to minimum phase but the filter length remains the same: Move all the zeros that are outside the unit circle to that they are radially inverted (and inside the unit circle). The order of the filter remains the same - so that doesn't meet your objective it appears. Here's a more interesting way: 1) Start by designing a linear phase FIR filter that has all-positive frequency response. This can be done as follows: a) Find a regular FIR filter with all stop-band ripples equal. b) Now find a new FIR filter with the stop-band specification equal to the peak ripple from (a). (This may be called a "pass band" with different "gain") c) If (b) doesn't have double zeros in the stop band, adjust the stop-band specification until the filter has double zeros in the stop band(s). All the while keeping the order of the filter constant..... [If you can get "inside" the design algorithm then you can change the stopband specification at each iteration of the Remez exchange algorithm by increasing the stopband specification by the minimum of the selected peaks. At first, the minimum of the selected error peaks is less than the maximum - so the stopband magnitude is increased such that it is inside the peak-to-peak range of the ripple. Eventually, the minimum of the peaks becomes equal to the maximum of the peaks, the negative peaks equal the positive peaks and the stop band error varies from zero to 2 times the peak error.] 2) Take the square root of this filter. Since the zeros show up in constellations of 4, an easy way to do that is to remove all of the zeros that are outside the unit circle and to make the zeros on the unit circle single zeros instead of the double zeros that were created above. This is a minimum phase filter with roughly half the coefficients. The only thing to note in this process is that the peak stopband error is the square root of the original. So, 1.0*e^-6 becomes 1.0*e^-3. This means you may have to design first for better stopband ripple in order to arrive at what you need. The method of defining an all-positive filter response and taking the square root is due to Schuessler. The method of modifying the stopband specification inside a Remez exchange algorithm is something I developed to enable working with Schuessler's process. I didn't ever prove convergence using this modified Remez exchange method - but it's never failed either! Fred
Reply by Fred Marshall March 1, 20042004-03-01
"v" <yes@i.am> wrote in message news:4043376A.8CD98C21@i.am...

.....................
> And here the algorithms like REMEZ > often can give very bad results > because it optimizes not the max. error > but the mean square error.
......................>
> A.Ser. >
A. Ser., Well, I suppose that one could label *any* algorithm improperly. However, the Remez exchange algorithm *is* directed at minimizing the maximum error. Now, I can't vouch for what may be in Matlab or any other "collection" with such labels attached. Nonetheless, I've never seen a counter example. Fred