DSPRelated.com
Forums

Halfband FIR design with alternating zeros

Started by I. R. Khan June 9, 2005
I tried matlab 'remez.m' and several online applets to design equiripple 
halfband lowpass FIR filters. All of them give the filter coeffiecients, 
which are all non-zero. I was expecting that every third coefficient in 
halfband designs should be zero. Is this correct? If yes, then is there some 
online free software that gives such results for equiripple filters. I have 
found such designs for maxflat filters but not for equiripple filters.

Thank you.

Ishtiaq. 

I. R. Khan wrote:
> I tried matlab 'remez.m' and several online applets to design equiripple > halfband lowpass FIR filters. All of them give the filter coeffiecients, > which are all non-zero. I was expecting that every third coefficient in > halfband designs should be zero. Is this correct? If yes, then is there > some online free software that gives such results for equiripple > filters. I have found such designs for maxflat filters but not for > equiripple filters. > > Thank you. > > Ishtiaq.
No it is not correct, every second (except the middle one) shold be zero. Matlab code for desinging halfband FIR filters is given below. -- Juha filterOrder=23; passbandEdge=0.4; [initH, ripple] = remez(filterOrder, [0 2*passbandEdge 1 1], ... [0.5 0.5 0 0], {256}) h((filterOrder+1):filterOrder:(2*filterOrder+1)) = [0.5 0]; h(1:2:end) = initH(1:end); freqz(h)
"I. R. Khan" <ir_khan@hotmail.com> wrote in message 
news:3gqbqaFdpiopU1@individual.net...
>I tried matlab 'remez.m' and several online applets to design equiripple >halfband lowpass FIR filters. All of them give the filter coeffiecients, >which are all non-zero. I was expecting that every third coefficient in >halfband designs should be zero. Is this correct? If yes, then is there >some online free software that gives such results for equiripple filters. I >have found such designs for maxflat filters but not for equiripple filters. > > Thank you. > > Ishtiaq.
There's a PC executable that I keep posted here: ftp://ftp.mission-systems-inc.com/outgoing/Halfband/ Fred
I. R. Khan wrote:
> I tried matlab 'remez.m' and several online applets to design equiripple > halfband lowpass FIR filters. All of them give the filter coeffiecients, > which are all non-zero. I was expecting that every third coefficient in > halfband designs should be zero. Is this correct? If yes, then is there > some online free software that gives such results for equiripple > filters. I have found such designs for maxflat filters but not for > equiripple filters.
Middle coefficient non zero, but all other odd-numbered coefficients. Round-off errors usually cause computer programs to return very small numbers instead of zero. Ignore those. Jerry -- Engineering is the art of making what you want from things you can get. &macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;
Fred Marshall wrote:
...
> There's a PC executable that I keep posted here: > > ftp://ftp.mission-systems-inc.com/outgoing/Halfband/
Fred, thanks! Andor
Thanks Juha for the code.

I don't understand some statements. Is [0 2*passbandEdge 1 1] in remez 
statement correct?
Could you please send the coefficients for one or two different small 
orders, say for filterOrder = 5 and/or 10, so that I can match it with the 
results I am getting.

For filterOrder = 6, I get
h = [0.0755, 0, -0.0433, 0, 0.0495, 0.5000,  0.4483, 0, 0.0495, 0, -0.0433, 
0,  0.0755}
and reponse does not seem to be a halfband lowpass.

Ishtiaq.

> filterOrder=23; passbandEdge=0.4; > [initH, ripple] = remez(filterOrder, [0 2*passbandEdge 1 1], ... > [0.5 0.5 0 0], {256}) > h((filterOrder+1):filterOrder:(2*filterOrder+1)) = [0.5 0]; > h(1:2:end) = initH(1:end); > freqz(h)
Fred, thanks for the code.

Just a comment. While matlab can easily import coefficients generated by the 
code, Mathematica has some problems due to use of E. You might consider 
using plane %g format for lower order designs.

Ishtiaq.


"Fred Marshall" <fmarshallx@remove_the_x.acm.org> wrote in message 
news:h9ydnfZxKcaFyTXfRVn-iQ@centurytel.net...

> > There's a PC executable that I keep posted here: > > ftp://ftp.mission-systems-inc.com/outgoing/Halfband/ > > Fred > > >
in article 5JqdnR5udeMrETXfRVn-vw@rcn.net, Jerry Avins at jya@ieee.org wrote
on 06/09/2005 14:38:

> Middle coefficient non zero, but all other odd-numbered coefficients. > Round-off errors usually cause computer programs to return very small > numbers instead of zero. Ignore those.
also, make awfully damn sure that the gain about Nyquist/2 is perfectly odd symmetry and that the weighting about Nyquist/2 is perfectly even symmetry. then you will get virtually half of the FIR coefs to be very small numbers that, as Jerry sez, should be set to zero. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
I. R. Khan wrote:
> Thanks Juha for the code. > > I don't understand some statements. Is [0 2*passbandEdge 1 1] in remez > statement correct? > Could you please send the coefficients for one or two different small > orders, say for filterOrder = 5 and/or 10, so that I can match it with > the results I am getting. > > For filterOrder = 6, I get > h = [0.0755, 0, -0.0433, 0, 0.0495, 0.5000, 0.4483, 0, 0.0495, 0, > -0.0433, 0, 0.0755} > and reponse does not seem to be a halfband lowpass. > > Ishtiaq. >
Jes, it is valid statement. The peculiar syntax is due to reason that there was a bug in matlab 6.5.1. Btw the "filterOrder" shoud be odd number. -- BR, Juha
>> ver
------------------------------------------------------------------------------------- MATLAB Version 6.5.1.199709 (R13) Service Pack 1 Operating System: OSF1 V5.1 732 alpha Java VM Version: Java 1.1.8-5 from Compaq Computer Corp. -------------------------------------------------------------------------------------
>> %%%%% Correct design in matlab 6.5.1 >> clear all >> filterOrder=5; passbandEdge=0.4; >> [initH, ripple] = remez(filterOrder, [0 2*passbandEdge 1 1], ...
[0.5 0.5 0 0], {256}) initH = 0.0537 -0.0915 0.3132 0.3132 -0.0915 0.0537 ripple = 0.0509
>> h((filterOrder+1):filterOrder:(2*filterOrder+1)) = [0.5 0]
h = Columns 1 through 7 0 0 0 0 0 0.5000 0 Columns 8 through 11 0 0 0 0
>> h(1:2:end) = initH(1:end)
h = Columns 1 through 7 0.0537 0 -0.0915 0 0.3132 0.5000 0.3132 Columns 8 through 11 0 -0.0915 0 0.0537
>> %%%%%% This doesn't work on 6.5.1 (an excellent error message!!) >> filterOrder=5; passbandEdge=0.4; >> [initH, ripple] = remez(filterOrder, [0 2*passbandEdge], [0.5 0.5],
{256}) at the Nyquist frequency. The order is being increased by one.ero
> In /commercial/matlabr13.1/toolbox/signal/signal/remez.m at line 204
*** FAILURE TO CONVERGE *** Probable cause is machine rounding error. Number of iterations = 3 Warning: If the number of iterations exceeds 3, the design may be correct, but should be verified with freqz. initH = 0.0000 0 0.0000 0.5000 0.0000 0 0.0000 ripple = 1.7270e-17
>> %%%%% But it works on matlab 7.0 >> ver
------------------------------------------------------------------------------------- MATLAB Version 7.0.4.352 (R14) Service Pack 2 Operating System: Linux 2.4.29-pre2-cs1 #2 SMP Mon Dec 20 11:59:34 EET 2004 i686 Java VM Version: Java is not enabled -------------------------------------------------------------------------------------
>> >> filterOrder=5; passbandEdge=0.4; >> [initH, ripple] = remez(filterOrder, [0 2*passbandEdge], [0.5 0.5],
{256}) initH = 0.0537 -0.0915 0.3132 0.3132 -0.0915 0.0537 ripple = 0.0509
>> h((filterOrder+1):filterOrder:(2*filterOrder+1)) = [0.5 0]
h = 0 0 0 0 0 0.5000 0 0 0 0 0
>> h(1:2:end) = initH(1:end)
h = 0.0537 0 -0.0915 0 0.3132 0.5000 0.3132 0 -0.0915 0 0.0537
"I. R. Khan" <ir_khan@hotmail.com> wrote in message 
news:3gsb31Fe2bocU1@individual.net...
> Thanks Juha for the code. > > I don't understand some statements. Is [0 2*passbandEdge 1 1] in remez > statement correct?
I found this statement interesting. It sugggests something I'd not tried using the Parks-McClellan program - which is to define a zero-length "band" with "x x", in this case "1 1". I tried these same parameters in the P-M program and it worked. The remez program in Matlab uses the normalization of fs/2=1 instead of the more common fs=1. Since the passbandEdge is stated as 0.4 (which is typical for halfband filters), meaning 0.4*fs, that's why it shows up as 2*passbandEdge. I also found the algorithmic approach interesting. First, the minimax filter response above is found for a filter of odd order / even length. Now, a filter of even length is either registered in time so that it's not even around zero or, to make it symmetrical around zero time one has to double the sample rate of the filter by inserting zeros. I like to do this so that the FFT is purely real making it easier to examine. So, this short bit of code first finds an antisymmetric lowpass filter of even length. Because the filter is even in time and of even length, the frequency response is only made of of odd harmonics of cosines. This means that the frequency response at fs/2 is zero (all odd cosines are zero at fs/2). Then, it intersperses N-1 zeros in time to double the sample rate as above. This results in frequency response zeros at the new fs/4 and 3fs/4 and a frequency response with zero integral. Then, it adds 0.5 which makes the antisymmetric response equal to 0.5 at fs/2 and 3fs/4, the passband value wiggling around 1.0 and the stopband value wiggling identically around zero(because the stopband is a mirror image of the passband). Probably this is the same trick for using P-M or "remez" rather directly for halfband filters that fred harris mentions in his book. The filter is symmetrical in time so length 24 gives 12 degrees of freedom and the passband ripple peaks count to N+1 or 13 as expected. Because the stated desired response at fs/2 is zero, the same as what the resulting response will be for a symmetrical, even-length filter, there is no error peak at fs/2. Fred