IIR filter implementation in C

Hi all, 

I have designed an algorithm in MATLAB, and am currently trying to
translate it into C code so that I can run it in real time. At the moment
I am stuck on the IIR filter. I am trying to implement the function 
[B,A] = cheby1(4,1,[80/48000,16000/48000],'bandpass')

I have found two promising sources online that provide source code for an
IIR filter in C:
1) http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html 
2) http://www.csois.usu.edu/people/yqchen/filter.c/F.C 

Both sources seem to work reasonably well for LPF cases (source two
provides a sample data set that shows how close the F.c function is to
filter.m), however, when trying to implement a BPF (or even a HPF), both
sources seem to die. 

For source 1, I have tried designing the filter using the website, and
then using the code that it outputs (with slight modification for
testability). Most of the time the resultant signal just oscillates
unstably. For source 2, I have designed different filters in MATLAB and
then copy/pasted the coefficients into the code (and adjusted the ORDER
variable accordingly). When using the filter coefficients generated from
cheby1(4,1,[80/48000,16000/48000],'bandpass'), the resultant signal grows
exponentially from 0 to 2e37. 

I have been unable to determine why either of these sources seems to
produce weird results for my BPF. I have made sure to use data sizes that
are powers of 2. I was under the impression that the recursion formula was
the same regardless of the filter coefficients, and so simply changing the
coefficients SHOULD be enough to change the behaviour of the filter code.
Apparently this isn�t the case. Can anyone see something in the code that
would need to be changed specifically for a BPF?? 

I have found some other code under GPL as well, however most of them make
use of the mkfilter routine (source 1), and so they have no worked
properly either. Ideally I would like to use source 2, due to its
compactness and simplicity. 

Any help that someone can offer would be more than greatly appreciated.

benr wrote:
> Hi all, 
> 
> I have designed an algorithm in MATLAB, and am currently trying to
> translate it into C code so that I can run it in real time. At the moment
> I am stuck on the IIR filter. I am trying to implement the function 
> [B,A] = cheby1(4,1,[80/48000,16000/48000],'bandpass')
> 
> I have found two promising sources online that provide source code for an
> IIR filter in C:
> 1) http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html 
> 2) http://www.csois.usu.edu/people/yqchen/filter.c/F.C 
> 
> Both sources seem to work reasonably well for LPF cases (source two
> provides a sample data set that shows how close the F.c function is to
> filter.m), however, when trying to implement a BPF (or even a HPF), both
> sources seem to die. 
> 
> For source 1, I have tried designing the filter using the website, and
> then using the code that it outputs (with slight modification for
> testability). Most of the time the resultant signal just oscillates
> unstably. For source 2, I have designed different filters in MATLAB and
> then copy/pasted the coefficients into the code (and adjusted the ORDER
> variable accordingly). When using the filter coefficients generated from
> cheby1(4,1,[80/48000,16000/48000],'bandpass'), the resultant signal grows
> exponentially from 0 to 2e37. 
> 
> I have been unable to determine why either of these sources seems to
> produce weird results for my BPF. I have made sure to use data sizes that
> are powers of 2. I was under the impression that the recursion formula was
> the same regardless of the filter coefficients, and so simply changing the
> coefficients SHOULD be enough to change the behaviour of the filter code.
> Apparently this isn�t the case. Can anyone see something in the code that
> would need to be changed specifically for a BPF?? 
> 
> I have found some other code under GPL as well, however most of them make
> use of the mkfilter routine (source 1), and so they have no worked
> properly either. Ideally I would like to use source 2, due to its
> compactness and simplicity. 
> 
> Any help that someone can offer would be more than greatly appreciated.

I'm a stranger to all your programs; discount my guess accordingly. I'll 
lay (small) odds that the order and interpretation of the coefficients 
are different from what you assume. In some developments, even the signs 
are reversed.

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

benr wrote:

> Hi all, 
> 
> I have designed an algorithm in MATLAB, and am currently trying to
> translate it into C code so that I can run it in real time. At the moment
> I am stuck on the IIR filter. I am trying to implement the function 
> [B,A] = cheby1(4,1,[80/48000,16000/48000],'bandpass')
> 
> I have found two promising sources online that provide source code for an
> IIR filter in C:
> 1) http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html 
> 2) http://www.csois.usu.edu/people/yqchen/filter.c/F.C 
> 
> Both sources seem to work reasonably well for LPF cases (source two
> provides a sample data set that shows how close the F.c function is to
> filter.m), however, when trying to implement a BPF (or even a HPF), both
> sources seem to die. 
> 
> For source 1, I have tried designing the filter using the website, and
> then using the code that it outputs (with slight modification for
> testability). Most of the time the resultant signal just oscillates
> unstably. For source 2, I have designed different filters in MATLAB and
> then copy/pasted the coefficients into the code (and adjusted the ORDER
> variable accordingly). When using the filter coefficients generated from
> cheby1(4,1,[80/48000,16000/48000],'bandpass'), the resultant signal grows
> exponentially from 0 to 2e37. 
> 
> I have been unable to determine why either of these sources seems to
> produce weird results for my BPF. I have made sure to use data sizes that
> are powers of 2. I was under the impression that the recursion formula was
> the same regardless of the filter coefficients, and so simply changing the
> coefficients SHOULD be enough to change the behaviour of the filter code.
> Apparently this isn�t the case. Can anyone see something in the code that
> would need to be changed specifically for a BPF?? 
> 
> I have found some other code under GPL as well, however most of them make
> use of the mkfilter routine (source 1), and so they have no worked
> properly either. Ideally I would like to use source 2, due to its
> compactness and simplicity. 
> 
> Any help that someone can offer would be more than greatly appreciated.
> 
> 
> 
What are the pole locations in your filters?  With direct-form 2nd-order 
filters you need data paths that are wide enough not only for the input 
data, but to adequately separate the pole positions from 1, twice.  So 
if the absolute value of the filter poles is 0.999 you need to have the 
data width + 20 bits (-2log2(1-0.999)).  If you're trying to implement a 
narrow bandpass this could be tripping you up.

-- 

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com

Posting from Google?  See http://cfaj.freeshell.org/google/

Thanx for replying so quickly guys. Turns out (after much experimentation)
that the sources work for order 1 and order 2 fitlers, but not for my
order 4 filter (even worse, the BPF is order 8). 

I guess it is something to do with the precision of the filter
coefficients, leading to unstable filter behaviour and oscillation. 

So I have solved the problem by implementing lower order filters. They are
still not completely accurate, but close enough for the moment. A little
more tweaking and I should be there. 

Thank you again

benr wrote:
> Thanx for replying so quickly guys. Turns out (after much experimentation)
> that the sources work for order 1 and order 2 fitlers, but not for my
> order 4 filter (even worse, the BPF is order 8). 
> 
> I guess it is something to do with the precision of the filter
> coefficients, leading to unstable filter behaviour and oscillation. 
> 
> So I have solved the problem by implementing lower order filters. They are
> still not completely accurate, but close enough for the moment. A little
> more tweaking and I should be there. 
> 
> Thank you again

A straight implementation of a high order IIR is almost always a 
disaster. To maintain any accuracy you need to split the filter into 
smaller chunks, and bi-quadratic sections are the sweet spot for doing 
that. Hence, most stuff you will see about IIR filters implements them 
as a chain of bi-quadratic sections.

There are 2 basic forms of bi-quad implementaion, and several variants, 
with much recent discussion about their merits. In fact, any time you 
come to this news group you will find there has been much recent 
discussion about their merits. :-)

Regards,
Steve

Perhaps you can save yourself aggravation by simply storing the
coefficients of your filter (B and A) as static data. Of course this
assumes that the characteristic of your filter does not need to change
during run time.