iir filter implementation

Started by Kuppala Sridhar March 19, 2007
hai,

iam trying to implement the iir filter on adsp21369 ezkit.the sampling rateis 48khz.
i have implemented the filter using second order sections
w[n]=x[n]+a1*w[n-1]+a2*w[n-2];
y[n]=w[n]+b1*w[n-1]+b2*w[n-2];

i have followed the example stereo cascade iir filter present in the analog devices website.but it is not working.
iam having doubt in the length of delayline.
please help me in implementing the iir filter.

---------------------------------
Heres a new way to find what you're looking for - Yahoo! Answers
On Mon, 19 Mar 2007, Kuppala Sridhar wrote:

> hai,
>
> iam trying to implement the iir filter on adsp21369 ezkit.the sampling rateis 48khz.
> i have implemented the filter using second order sections
> w[n]=x[n]+a1*w[n-1]+a2*w[n-2];
> y[n]=w[n]+b1*w[n-1]+b2*w[n-2];
>
> i have followed the example stereo cascade iir filter present in the analog devices website.but it is not working.
> iam having doubt in the length of delayline.
> please help me in implementing the iir filter.

Howdy Kuppala,

What values do you have for a1, a2, b1 and b2? As a first attempt to test
your filters, just try one biquad at a time. Also put in known
frequencies - use a table look up to fill a buffer and just cycle thru it.
Each biquad should reduce the amplitude according to where you are in the
pass band of that filter. I have found that just having one number wrong
is all it takes to totally screw up a filter, and finding it can be a
pain. But if you do one biquad at a time, it pops out pretty fast.

Patience, persistence, truth,
Dr. mike
hai mike,
thanks for u r reply.
i have also implemented with single biquad but this is not responding the way the frquency response in the matlab was.
how to decide a1,a2,b1,b2 from the sos matrix generated in matlab.we have taken these values from the directform II canonical form realized model generated in matlab.some times even though in fda tool it is showing that the filter is stable but it is behaving unstable(dline buf filled with all 0xFFFFFFFF)

Mike Rosing wrote:
On Mon, 19 Mar 2007, Kuppala Sridhar wrote:

> hai,
>
> iam trying to implement the iir filter on adsp21369 ezkit.the sampling rateis 48khz.
> i have implemented the filter using second order sections
> w[n]=x[n]+a1*w[n-1]+a2*w[n-2];
> y[n]=w[n]+b1*w[n-1]+b2*w[n-2];
>
> i have followed the example stereo cascade iir filter present in the analog devices website.but it is not working.
> iam having doubt in the length of delayline.
> please help me in implementing the iir filter.

Howdy Kuppala,

What values do you have for a1, a2, b1 and b2? As a first attempt to test
your filters, just try one biquad at a time. Also put in known
frequencies - use a table look up to fill a buffer and just cycle thru it.
Each biquad should reduce the amplitude according to where you are in the
pass band of that filter. I have found that just having one number wrong
is all it takes to totally screw up a filter, and finding it can be a
pain. But if you do one biquad at a time, it pops out pretty fast.

Patience, persistence, truth,
Dr. mike
On Wed, 21 Mar 2007, Kuppala Sridhar wrote:

> hai mike,
> thanks for u r reply.
> i have also implemented with single biquad but this is not responding the way the frquency response in the matlab was.
> how to decide a1,a2,b1,b2 from the sos matrix generated in matlab.we have taken these values from the directform II canonical form realized model generated in matlab.some times even though in fda tool it is showing that the filter is stable but it is behaving unstable(dline buf filled with all 0xFFFFFFFF)
>
> Mike Rosing wrote:
> On Mon, 19 Mar 2007, Kuppala Sridhar wrote:
>
> > hai,
> >
> > iam trying to implement the iir filter on adsp21369 ezkit.the sampling rateis 48khz.
> > i have implemented the filter using second order sections
> > w[n]=x[n]+a1*w[n-1]+a2*w[n-2];
> > y[n]=w[n]+b1*w[n-1]+b2*w[n-2];

Are you using floating point or fixed point? If fixed point, check the
scaling. Matlab doesn't know about your processor, you may have to do a
divide by 2^16 or 2^32 after each multiply-accumulate.
hai mike,
> thanks for u r reply.
> i have also implemented with single biquad but this is not responding the way the frquency response in the matlab was.

Hello Kuppala,

here are some additional hints, which might be helpful or not:

1) be aware that Matlab specifies the coefficient differently from how it is usually done in books or code libraries. Usually they differ in sign.
Check the algorithm thoroughly.

If you select coefficients with easy numbers like 0.0 0.5 or 1.0 where appropriate, you can calculate the biquad manually and check if you get the same results.

2) Depending on the cut-off frequency and the sampling frequency, you might run into wrongly calculated coefficients (precision issue).
For one biquad and a low-pass filter stage, keeping the relation
Fc/Fs > 0.001 should be safe.
Check the coefficients. Be cautious, if the values are very close to 0.0 or 1.0 or -1.0

3) Check what your compiler does. If it calculates the biquad based on 16 bit wide values, the algorithm will fail or at least provide a distorted result. Even 32bit wide values might be unsufficient, depending on type and algorithm.
Keep in mind, that IIR filters do infinite recursions on the calculated values. This means, that errors might increase at every cycle.
If you can watch the algorithm in a debugger, look for sign changes (values around zero) and maybe saturation or rollover effects at the border of the value range (>=1.0, <=-1,0 ).

Good luck.
Bernhard
Hello guys,

I'd like to emphasize the first remark of Bernhard, specially the sign of the recursive coefficientes (a_1, a_2). I've often found that simply changing their sign makes it work, because of the structure that each routine (Matlab's or DSP's) assume.

Regards,

Jaime Andr Aranguren Cardona
j...@ieee.org
j...@computer.org

----- Original Message ----
From: "H...@Foerstergroup.de"
To: a...
Sent: Thursday, March 22, 2007 1:30:28 AM
Subject: [adsp] Re: iir filter implementation

hai mike,
> thanks for u r reply.
> i have also implemented with single biquad but this is not responding the way the frquency response in the matlab was.

Hello Kuppala,

here are some additional hints, which might be helpful or not:

1) be aware that Matlab specifies the coefficient differently from how it is usually done in books or code libraries. Usually they differ in sign.
Check the algorithm thoroughly.

If you select coefficients with easy numbers like 0.0 0.5 or 1.0 where appropriate, you can calculate the biquad manually and check if you get the same results.

2) Depending on the cut-off frequency and the sampling frequency, you might run into wrongly calculated coefficients (precision issue).
For one biquad and a low-pass filter stage, keeping the relation
Fc/Fs > 0.001 should be safe.
Check the coefficients. Be cautious, if the values are very close to 0.0 or 1.0 or -1.0

3) Check what your compiler does. If it calculates the biquad based on 16 bit wide values, the algorithm will fail or at least provide a distorted result. Even 32bit wide values might be unsufficient, depending on type and algorithm.
Keep in mind, that IIR filters do infinite recursions on the calculated values. This means, that errors might increase at every cycle.
If you can watch the algorithm in a debugger, look for sign changes (values around zero) and maybe saturation or rollover effects at the border of the value range (>=1.0, <=-1,0 ).

Good luck.
Bernhard

____________________________________________________________________________________
Looking for earth-friendly autos?
Browse Top Cars by "Green Rating" at Yahoo! Autos' Green Center.
http://autos.yahoo.com/green_center/
hello Bernhard,
as per ur suggestions, i have made changes with respect to the sign change of the filter coefficients.But even then my frequency response is not proper.please see the code below whether any bug is there.

#define FRAMESIZE 1 /* stereo */
#define SECTIONS 1 /* Number of second-order sections (biquads) */
/* Each channel will be filtered by SECTIONS biquads */
/***********************************************************************/
.section /dm seg_dmda;

.var inbuf[FRAMESIZE]; /* All of the channels in each frame are processed */
/* when the last channel arrives and the */

.var outbuf[FRAMESIZE]; /* double-buffering feeds the SPORT as needed */
/* tag data to check if ever assigned a value */

.var delaybuf[SECTIONS*4]; /* Holds intermediate variables (state) for the */
/* IIR. Each of the two compute units needs */
/* storage for w[x-1] and w[x-2] for each */
/* biquad section. */

/***********************************************************************/
.section /pm seg_pmda;

/* Coefficients are ordered a2,a1,b2,b1 and coefficients for PEx are */
/* immediately followed by those for PEy. Repeat the sequence for each */
/* successive cascaded biquad pair. */

.var coefficients[SECTIONS*8] {
1.341081473676422 , - 0.571654398453466,
1.994354460540368 , 0.999999999999998 ,
};

/***********************************************************************/
.section /pm seg_pmco;
_initializeAlgorithmDataStructures:

m1=1; /* modifiers=2 for SIMD reads & writes */
m2=1;
m8=1;

/* Not using circular addressing with any of the pointers, so */
/* we must ensure that it is disabled by setting length = 0 */

l3=0; /* i3 points to inbuf */
l4=0; /* i4 points to outbuf */

/* i0 and i1 are pointer into delaybuf */
/* and track w[n-1] and w[n-2] */

l0=4; /* w[n-1], w[n-2] for PEx, PEy */
l1=0;

l8=0; /* i8 points to the coefficient buffer */

/* zero the delay line */

r0=SECTIONS;
b0aybuf;
f2=0;

lcntr=r0, do clear until lce; /* for each section, do: */
dm(i0,m2) /* clear w[n-1] and w[n-2] */
clear: dm(i0,m2)

rts;

_initializeAlgorithmDataStructures.end:
cascaded_biquad:
bit set mode1 CBUFEN;// // Enable SIMD mode
b0 = delaybuf;
b1 = b0;
b3 = inbuf;
b4 = outbuf;
b9 = coefficients;
r0 = SECTIONS;

f8=dm(i3,m1); // read inbuf
r12=r12 xor r12, f2=dm(i0,m1), f4=pm(i8,m8);
lcntr=r0, do quads until lce;
f12f4, f8f12, f3=dm(i0,m1), f4=pm(i8,m8);
f12f4, f8f12, dm(i1,m1) f4=pm(i8,m8);
f12f4, f8f12, f2=dm(i0,m1), f4=pm(i8,m8);
quads: f12f4, f8f12, dm(i1,m1) f4=pm(i8,m8);
f8f12;
rts (db);
dm(i4,m1)
bit clr mode1 CBUFEN ;// // disable SIMD mode

_Cascaded_IIR_Filter_SIMD.end:

H...@Foerstergroup.de wrote: hai mike,
> thanks for u r reply.
> i have also implemented with single biquad but this is not responding the way the frquency response in the matlab was.

Hello Kuppala,

here are some additional hints, which might be helpful or not:

1) be aware that Matlab specifies the coefficient differently from how it is usually done in books or code libraries. Usually they differ in sign.
Check the algorithm thoroughly.

If you select coefficients with easy numbers like 0.0 0.5 or 1.0 where appropriate, you can calculate the biquad manually and check if you get the same results.

2) Depending on the cut-off frequency and the sampling frequency, you might run into wrongly calculated coefficients (precision issue).
For one biquad and a low-pass filter stage, keeping the relation
Fc/Fs > 0.001 should be safe.
Check the coefficients. Be cautious, if the values are very close to 0.0 or 1.0 or -1.0

3) Check what your compiler does. If it calculates the biquad based on 16 bit wide values, the algorithm will fail or at least provide a distorted result. Even 32bit wide values might be unsufficient, depending on type and algorithm.
Keep in mind, that IIR filters do infinite recursions on the calculated values. This means, that errors might increase at every cycle.
If you can watch the algorithm in a debugger, look for sign changes (values around zero) and maybe saturation or rollover effects at the border of the value range (>=1.0, <=-1,0 ).

Good luck.
Bernhard

---------------------------------
Heres a new way to find what you're looking for - Yahoo! Answers
Hi Kuppala,

I didn't check all of your code.
However, I found one bug, where you try to enable SIMD mode:

> bit set mode1 CBUFEN;// // Enable SIMD mode
> bit clr mode1 CBUFEN ;// // disable SIMD mode

I would switch this on/off with the bit PEYEN in the same register.

If this doesn't solve your issue, use simpler values for the coefficients
and step through the biquad, watching the values how they change.
For a simple lowpass filter, you should be able to prepare this task
on a sheet of paper or in a spreadsheet program.
Then you might want to use the load memory function, thus restarting every run with the same preconditions.
Always keep in mind, that it's very difficult to code something which you don't fully understand, because there might be coding issues, design flaws and lots of other aspects causing head-ache.
So make sure that you really understand what every line of code should produce, then check thoroughly if it does.
Step by step.

Hint: I'd propose that you don't enable SIMD mode unless all of your code
is already behaving well in SISD mode.
(At least for newcomers) it's quite tricky what might happen once SIMD is enabled.
Bernhard
Dear friends,

Kindly let me know whether it is possible to process 2-D arrays directly ( I mean without converting the data to 1-D variables) in visual DSP directly.

Also kindly let me know how to go about plotting 2-d data in visual dsp.

eagerly waiting for your reply.

regards.

amit

________________________________
From: Kuppala Sridhar
To: a...
Sent: Saturday, 24 March 2007 2:19 AM
Subject: Re: [adsp] Re: iir filter implementation

hello Bernhard,
as per ur  suggestions,  i have made changes with respect to the sign change of the filter coefficients. But even then my frequency response is not proper.please see the code below whether any bug is there.

#define FRAMESIZE 1   /* stereo  */
#define SECTIONS 1    /* Number of second-order sections (biquads)         */
                      /* Each channel will be filtered by SECTIONS biquads */
/*********** ********* ********* ********* ********* ********* ********* ******/
.section /dm seg_dmda;

.var  inbuf[FRAMESIZE] ;  /* All of the channels in each frame are processed */
                         /*   when
the last channel arrives and the         */

.var  outbuf[FRAMESIZE] ; /*   double-buffering feeds the SPORT as needed    */
                         /*  tag data to check if ever assigned a value     */

.var  delaybuf[SECTIONS* 4]; /* Holds intermediate variables (state) for the */
                            /*   IIR.  Each of the two compute units needs  */
                            /*   storage for w[x-1] and w[x-2] for
each     */
                            /*   biquad section.                            */

/*********** ********* ********* ********* ********* ********* ********* ******/
.section /pm seg_pmda;

 /* Coefficients are ordered a2,a1,b2,b1 and coefficients for PEx are      */
 /*   immediately followed by those for PEy.  Repeat the sequence for each */
 /*   successive cascaded biquad
pair.                                     */

.var   coefficients[ SECTIONS* 8] {
     1.341081473676422 ,  - 0.571654398453466,
1.994354460540368 ,   0.999999999999998 ,     
};

/*********** ********* ********* ********* ********* ********* ********* ******/
.section /pm seg_pmco;
_initializeAlgorith mDataStructures:

        m1=1;  /* modifiers=2 for SIMD reads & writes */
        m2=1;
        m8=1;

    /* Not using circular addressing with any of the pointers, so */
    /*  we must ensure that it is disabled by setting length = 0 
*/

        l3=0;  /* i3 points to inbuf */
        l4=0;  /* i4 points to outbuf */

        /* i0 and i1 are pointer into delaybuf */
        /*    and track w[n-1] and w[n-2]      */

        l0=4;  /* w[n-1], w[n-2] for PEx, PEy */
        l1=0;

        l8=0;  /* i8 points to the coefficient buffer */

        /* zero the delay line */

        r0=SECTIONS;
        b0aybuf;
        f2=0;

        lcntr=r0, do clear until
lce;     /* for each section, do:       */
            dm(i0,m2)                  /*   clear w[n-1] and w[n-2]  */
clear:      dm(i0,m2)

    rts;

_initializeAlgorith mDataStructures. end:
cascaded_biquad:
    bit set mode1 CBUFEN;//           // Enable SIMD mode
    b0 = delaybuf;
    b1 = b0;
    b3 = inbuf;
    b4 = outbuf;
    b9 = coefficients;
    r0 = SECTIONS;

   
f8=dm(i3,m1) ;                               // read inbuf
    r12=r12 xor r12,    f2=dm(i0,m1) ,  f4=pm(i8,m8) ;
       lcntr=r0, do quads until lce;
            f12f4,         f8f12,       f3=dm(i0,m1) ,     f4=pm(i8,m8) ;
            f12f4,         f8f12,      dm(i1,m1)     f4=pm(i8,m8) ;
            f12f4,         f8f12,       f2=dm(i0,m1) ,    
f4=pm(i8,m8) ;
quads:      f12f4,         f8f12,       dm(i1,m1)     f4=pm(i8,m8) ;
    f8f12;
    rts (db);
    dm(i4,m1)
    bit clr mode1 CBUFEN ;//                  // disable SIMD mode

_Cascaded_IIR_ Filter_SIMD. end:

Holzmayer.Bernhard@ Foerstergroup. de wrote:
hai mike,
>> thanks for u
r reply.
>> i have also implemented with single biquad but this is not responding the way the frquency response in the matlab was.
>
>Hello Kuppala,
>
>here are some additional hints, which might be helpful or not:
>
>1) be aware that Matlab specifies the coefficient differently from how it is usually done in books or code libraries. Usually they differ in sign.
>Check the algorithm thoroughly.
>
>If you select coefficients with easy numbers like 0.0 0.5 or 1.0 where appropriate, you can calculate the biquad manually and check if you get the same results.
>
>2) Depending on the cut-off frequency and the sampling frequency, you might run into wrongly calculated coefficients (precision issue).
>For one biquad and a low-pass filter stage, keeping the relation
>Fc/Fs > 0.001 should be safe.
>Check the coefficients. Be cautious, if the values are very close to 0.0 or 1.0 or -1.0
>
>3) Check what your compiler does. If
it calculates the biquad based on 16 bit wide values, the algorithm will fail or at least provide a distorted result. Even 32bit wide values might be unsufficient, depending on type and algorithm.
>Keep in mind, that IIR filters do infinite recursions on the calculated values. This means, that errors might increase at every cycle.
>If you can watch the algorithm in a debugger, look for sign changes (values around zero) and maybe saturation or rollover effects at the border of the value range (>=1.0, <=-1,0 ).
>
>Good luck.
>Bernhard

________________________________
Here’s a new way to find what you're looking for - Yahoo! Answers
You can use the view->debug windows->image viewer to see the data from a
2d array directly. You can only process one memory location at a time,
but you can do a 2D array easily.

Patience, persistence, truth,
Dr. mike

On Thu, 26 Apr 2012, Amit Kumar Verma wrote:

> Dear friends,
>
> Kindly let me know whether it is possible to process 2-D arrays directly ( I mean without converting the data to 1-D variables) in visual DSP directly.
>
> Also kindly let me know how to go about plotting 2-d data in visual dsp.
>
> eagerly waiting for your reply.
>
> regards.
>
> amit
>
> ________________________________
> From: Kuppala Sridhar
> To: a...
> Sent: Saturday, 24 March 2007 2:19 AM
> Subject: Re: [adsp] Re: iir filter implementation
> hello Bernhard,
> as per ur  suggestions,  i have made changes with respect to the sign change of the filter coefficients. But even then my frequency response is not proper.please see the code below whether any bug is there.
>
> #define FRAMESIZE 1   /* stereo  */
> #define SECTIONS 1    /* Number of second-order sections (biquads)         */
>                       /* Each channel will be filtered by SECTIONS biquads */
> /*********** ********* ********* ********* ********* ********* ********* ******/
> .section /dm seg_dmda;
>
> .var  inbuf[FRAMESIZE] ;  /* All of the channels in each frame are processed */
>                          /*   when
> the last channel arrives and the         */
>
> .var  outbuf[FRAMESIZE] ; /*   double-buffering feeds the SPORT as needed    */
>                          /*  tag data to check if ever assigned a value     */
>
> .var  delaybuf[SECTIONS* 4]; /* Holds intermediate variables (state) for the */
>                             /*   IIR.  Each of the two compute units needs  */
>                             /*   storage for w[x-1] and w[x-2] for
> each     */
>                             /*   biquad section.                            */
>
> /*********** ********* ********* ********* ********* ********* ********* ******/
> .section /pm seg_pmda;
>
>  /* Coefficients are ordered a2,a1,b2,b1 and coefficients for PEx are      */
>  /*   immediately followed by those for PEy.  Repeat the sequence for each */
>  /*   successive cascaded biquad
> pair.                                     */
>
> .var   coefficients[ SECTIONS* 8] > {
>      1.341081473676422 ,  - 0.571654398453466,
> 1.994354460540368 ,   0.999999999999998 ,     
> };
>
> /*********** ********* ********* ********* ********* ********* ********* ******/
> .section /pm seg_pmco;
> _initializeAlgorith mDataStructures:
>
>         m1=1;  /* modifiers=2 for SIMD reads & writes */
>         m2=1;
>         m8=1;
>
>     /* Not using circular addressing with any of the pointers, so */
>     /*  we must ensure that it is disabled by setting length = 0 
> */
>
>         l3=0;  /* i3 points to inbuf */
>         l4=0;  /* i4 points to outbuf */
>
>         /* i0 and i1 are pointer into delaybuf */
>         /*    and track w[n-1] and w[n-2]      */
>
>         l0=4;  /* w[n-1], w[n-2] for PEx, PEy */
>         l1=0;
>
>         l8=0;  /* i8 points to the coefficient buffer */
>
>         /* zero the delay line */
>
>         r0=SECTIONS;
>         b0aybuf;
>         f2=0;
>
>         lcntr=r0, do clear until
> lce;     /* for each section, do:       */
>             dm(i0,m2)                  /*   clear w[n-1] and w[n-2]  */
> clear:      dm(i0,m2)
>
>     rts;
>
> _initializeAlgorith mDataStructures. end:
> cascaded_biquad:
>     bit set mode1 CBUFEN;//           // Enable SIMD mode
>     b0 = delaybuf;
>     b1 = b0;
>     b3 = inbuf;
>     b4 = outbuf;
>     b9 = coefficients;
>     r0 = SECTIONS;
>
>    
> f8=dm(i3,m1) ;                               // read inbuf
>     r12=r12 xor r12,    f2=dm(i0,m1) ,  f4=pm(i8,m8) ;
>        lcntr=r0, do quads until lce;
>             f12f4,         f8f12,       f3=dm(i0,m1) ,     f4=pm(i8,m8) ;
>             f12f4,         f8f12,      dm(i1,m1)     f4=pm(i8,m8) ;
>             f12f4,         f8f12,       f2=dm(i0,m1) ,    
> f4=pm(i8,m8) ;
> quads:      f12f4,         f8f12,       dm(i1,m1)     f4=pm(i8,m8) ;
>     f8f12;
>     rts (db);
>     dm(i4,m1)
>     bit clr mode1 CBUFEN ;//                  // disable SIMD mode
>
> _Cascaded_IIR_ Filter_SIMD. end:
>
> Holzmayer.Bernhard@ Foerstergroup. de wrote:
> hai mike,
>>> thanks for u
> r reply.
>>> i have also implemented with single biquad but this is not responding the way the frquency response in the matlab was.
>>
>> Hello Kuppala,
>>
>> here are some additional hints, which might be helpful or not:
>>
>> 1) be aware that Matlab specifies the coefficient differently from how it is usually done in books or code libraries. Usually they differ in sign.
>> Check the algorithm thoroughly.
>>
>> If you select coefficients with easy numbers like 0.0 0.5 or 1.0 where appropriate, you can calculate the biquad manually and check if you get the same results.
>>
>> 2) Depending on the cut-off frequency and the sampling frequency, you might run into wrongly calculated coefficients (precision issue).
>> For one biquad and a low-pass filter stage, keeping the relation
>> Fc/Fs > 0.001 should be safe.
>> Check the coefficients. Be cautious, if the values are very close to 0.0 or 1.0 or -1.0
>>
>> 3) Check what your compiler does. If
> it calculates the biquad based on 16 bit wide values, the algorithm will fail or at least provide a distorted result. Even 32bit wide values might be unsufficient, depending on type and algorithm.
>> Keep in mind, that IIR filters do infinite recursions on the calculated values. This means, that errors might increase at every cycle.
>> If you can watch the algorithm in a debugger, look for sign changes (values around zero) and maybe saturation or rollover effects at the border of the value range (>=1.0, <=-1,0 ).
>>
>> Good luck.
>> Bernhard
>>
>> ________________________________
> Here’s a new way to find what you're looking for - Yahoo! Answers