Digital IIR Antialias filters

Technical questions regarding the XTC tools and programming with XMOS.
User avatar
lilltroll
XCore Expert
Posts: 956
Joined: Fri Dec 11, 2009 3:53 am
Location: Sweden, Eskilstuna

Digital IIR Antialias filters

Post by lilltroll »

Test with IIR antialias filters here. For an example to be used in a multirate environment.
This test is when you need a very steep transitionband.
First out is a 10 order butterworth transformed to Second Order Sections - just to check how much roundoff noise we get with
int64 + int32*int32 => int64 (compared to double prec. without internal casting)

Code: Select all

#include <xs1.h>
#define BANKS 5

int X[BANKS][3]={{0,0,0},{0,0,0},{0,0,0},{0,0,0},{0,0,0}};
int Y[BANKS][2]={{0,0},{0,0},{0,0},{0,0},{0,0}};


int IIRfix64(int x){   
	const int B[BANKS][3]={ {8409818, 16819636, 8409818} , {8136995 ...
	const int A[BANKS][2]={ {1949207645, -1784570915} , {18859734 ...
	int Xin,h;
	unsigned int l,k;
	
	Xin=x<<6;
   for(k=0;k<BANKS;k++){
    h=0;l=0;
   {h, l} = macs(B[k][0], Xin , h, l);
   {h, l} = macs(B[k][1], X[k][1], h, l);
   {h, l} = macs(B[k][2], X[k][2], h, l);
   {h, l} = macs(A[k][0], Y[k][0]<<1, h, l);
   {h, l} = macs(A[k][1], Y[k][1], h, l);
   X[k][2]=X[k][1];
   X[k][1]=Xin;
   Y[k][1]=Y[k][0];
   Xin=h<<1;
   Y[k][0]=Xin;
   }
        
   return (Xin>>6);
}
The code is just for benchmarking the amount of roundoff noise - and is not optimized for computational speed.

The SOS sections looks like this with a -3 dB at 1 kHz @48 kHz fs
Butt10.png
Butt10.png (11.67 KiB) Viewed 4185 times
Butt10.png
Butt10.png (11.67 KiB) Viewed 4185 times
And the cumulative sections forms to this one.
Butt10cum.png
Butt10cum.png (12.47 KiB) Viewed 4185 times
Butt10cum.png
Butt10cum.png (12.47 KiB) Viewed 4185 times
One important thing here is that each cumulation never gets above 0 dB. Othervise clipping could occur in between the links. Scaling for fixed represenation can be done in many different ways. The one above is not the optimal one regarding to Hi Q-value polepars.
Last edited by lilltroll on Sun Mar 07, 2010 12:34 pm, edited 1 time in total.


Probably not the most confused programmer anymore on the XCORE forum.
User avatar
lilltroll
XCore Expert
Posts: 956
Joined: Fri Dec 11, 2009 3:53 am
Location: Sweden, Eskilstuna

Post by lilltroll »

The input signal to the filter was quantized to 24-bits resolution.
The total AC-roundoff noise in the filter is less than any common Audio-CODEC can achive.

Next up is the Elliptic filter - thats a more tricky one.
Attachments
Butt10_noise.png
Butt10_noise.png (6.37 KiB) Viewed 4182 times
Butt10_noise.png
Butt10_noise.png (6.37 KiB) Viewed 4182 times
Probably not the most confused programmer anymore on the XCORE forum.
User avatar
lilltroll
XCore Expert
Posts: 956
Joined: Fri Dec 11, 2009 3:53 am
Location: Sweden, Eskilstuna

Post by lilltroll »

For the elliptic filter much more care has been taken regarding the ordering of the 2nd order sections and the scaling. This filter will create a larger roundoff noise, but is also the steepest one possible for an given filter order.
Attachments
Elipp10cum.png
Elipp10cum.png (10.96 KiB) Viewed 4178 times
Elipp10cum.png
Elipp10cum.png (10.96 KiB) Viewed 4178 times
Probably not the most confused programmer anymore on the XCORE forum.
User avatar
lilltroll
XCore Expert
Posts: 956
Joined: Fri Dec 11, 2009 3:53 am
Location: Sweden, Eskilstuna

Post by lilltroll »

The error compared to double prec. - but the filter coeff is quantized to 32 bits prec. prior filtering with double.
(The scale of the error is 20 dB less compared to the Butterworth case)
Ellip10_noise.png
Ellip10_noise.png (7.36 KiB) Viewed 4174 times
Ellip10_noise.png
Ellip10_noise.png (7.36 KiB) Viewed 4174 times
Here it's worth to take a look at the spectral components of the error. Since the elliptic filter has a pole very close to 1 kHz - I expect an error maxima around 1 kHz.
Ellip10_PSD.png
Ellip10_PSD.png (8.77 KiB) Viewed 4173 times
Ellip10_PSD.png
Ellip10_PSD.png (8.77 KiB) Viewed 4173 times
Probably not the most confused programmer anymore on the XCORE forum.
User avatar
lilltroll
XCore Expert
Posts: 956
Joined: Fri Dec 11, 2009 3:53 am
Location: Sweden, Eskilstuna

Post by lilltroll »

Is it worth to add one more bit of prec. between each filter section on the cost of some extra instructions ?

Thus, change from

Code: Select all

Xin=h<<1

to

Code: Select all

Xin=(h<<1) + (l>>31);
Since I belive that this is the largest contributor to error since the accumulator is 64 bits long - lets compute the result.

The other question is - do I have more than 6 dB headroom in the accumulator even for worst case senario ?
Well - since I use the Sign mac I get one bit for free of headroom
int32*int32 = sign * uint31 * sign uint31 = sign * uint62

On the other hand I must be able to accumulate coef. B0-B2 without overflow.
Also I need headroom for the one bit Left Shift of Y(n) correcponding to coef. A1
Also, some of the filtersections sums to a little above 0 dB to avoid coef. overflow in the last filtersection with the nasty ploe-par. :geek:

:? Actually I do not believe that I can change from <<6 to <<7 for 24 bits input data.
But, lets make an trial and error

Some computing later...
<<7 is not possible - the filter overflows.
Back to <<6

... :idea: the overflow probably happend in the A1 mac. That could maybe be avoided with one extra instruction and some rewriting of the filter.

What about this instead :idea:

Code: Select all

	Xin=x<<6;
   for(k=0;k<BANKS;k++){
    h=0;l=0;
   {h, l} = macs(B[k][0], Xin , h, l);
   {h, l} = macs(B[k][1], X[k][1], h, l);
   {h, l} = macs(B[k][2], X[k][2], h, l);
   {h, l} = macs(A[k][0], Y[k][0], h, l);
   {h, l} = macs(A[k][1], Y[k][1]>>1, h, l);
   X[k][2]=X[k][1];
   X[k][1]=Xin;
   Y[k][1]=Y[k][0];
   Y[k][0]=(h<<2) + (l>>30);
   Xin=Y[k][0]>>1;
   }
Probably not the most confused programmer anymore on the XCORE forum.
User avatar
lilltroll
XCore Expert
Posts: 956
Joined: Fri Dec 11, 2009 3:53 am
Location: Sweden, Eskilstuna

Post by lilltroll »

Nice - Reduced the DC offset as well with a more correct rounding method of the LSB.
Ellip10_noise2.png
Ellip10_noise2.png (5.91 KiB) Viewed 4163 times
Ellip10_noise2.png
Ellip10_noise2.png (5.91 KiB) Viewed 4163 times
Lets check the diffence in the PSD.
Ellip10_PSD2.png
Ellip10_PSD2.png (6.6 KiB) Viewed 4163 times
Ellip10_PSD2.png
Ellip10_PSD2.png (6.6 KiB) Viewed 4163 times
When the signal energy is in the stopband - the noise gets even smaller.
Here a PSD over the last 3 seconds of the 10 s long chirp.
Ellip10_PSD710.png
Ellip10_PSD710.png (7.05 KiB) Viewed 4163 times
Ellip10_PSD710.png
Ellip10_PSD710.png (7.05 KiB) Viewed 4163 times
Probably not the most confused programmer anymore on the XCORE forum.
User avatar
Andy
Respected Member
Posts: 279
Joined: Fri Dec 11, 2009 1:34 pm

Post by Andy »

Interesting work lilltroll - do you plan to release these as part of your library?
User avatar
lilltroll
XCore Expert
Posts: 956
Joined: Fri Dec 11, 2009 3:53 am
Location: Sweden, Eskilstuna

Post by lilltroll »

Yes - I will, and I'm getting better on ASM now. Why ?

Well - I am using XTA in debug mode - but the debug mode was set to Optimization 0 (O0) :evil: by me then I played around with the mouse in XDE along time ago.

So when I checked the code in XTA I often realised that I could write the code better myself.

So I have had some ASM traning by looking at the compiled output in O0 and trying to write it better. Very good traning - since I hade something to look at - read the instruction list and learn.
Also I have learned more about how to write things in C - to not achieve stupid solutions at level O0. For an example I know that the way I have written the test code is stupid from an performance perspective. With X[] and Y[] outside the function, (instead of the use of references or pointers) - all values will always travel out to RAM and back to the registers between each arithmetic instruction.

But I will wait to write it smarter until I'm satiesfied with the numerical results. Next up is to studie the effects with dithering.

It's much harder to beat O3 ;)
Probably not the most confused programmer anymore on the XCORE forum.
User avatar
lilltroll
XCore Expert
Posts: 956
Joined: Fri Dec 11, 2009 3:53 am
Location: Sweden, Eskilstuna

Post by lilltroll »

An important question now - is there a significant linear correlation between the error and the signal e.g. Is the error mosty arising from a deviation in the frequency response due to limited coeff. prec. ? Such type of error is often a lesser problem compared to non linear errors - for an example in an antialias filter. The elliptic filter will ripple in the passband - in this example it ripples 0.1 dB (it's a design parameter) - if it instead ripples 0.1001 dB due to limited prec. that doesn't matter at all and that shouldn't be added as a roundoff problem in the analasis.

One way to check this out is to draw a spectrogram over the error.
If (only) one exponetial function can be seen in the plot - it will be the fundamental.
In this example the input signal was dithered to avoid higher partials.
Ellip_specgram.png
Ellip_specgram.png (55.9 KiB) Viewed 4154 times
Ellip_specgram.png
Ellip_specgram.png (55.9 KiB) Viewed 4154 times
It's hiding a very strong partial in the error :!:

If we take out a time slice at 5 s - 5.34 s and observe the PSD, it's easier to see the dB scale instead of tracking colours in the spectogram:
Ellip_psd5.png
Ellip_psd5.png (9.87 KiB) Viewed 4153 times
Ellip_psd5.png
Ellip_psd5.png (9.87 KiB) Viewed 4153 times
Maybe you are questioning why the roundoff noise has the shape it has?
Check out the last 2 filter sections in the filter!
ellip_sect.png
ellip_sect.png (9.38 KiB) Viewed 4152 times
ellip_sect.png
ellip_sect.png (9.38 KiB) Viewed 4152 times
Probably not the most confused programmer anymore on the XCORE forum.
User avatar
lilltroll
XCore Expert
Posts: 956
Joined: Fri Dec 11, 2009 3:53 am
Location: Sweden, Eskilstuna

Post by lilltroll »

Next up is the antialias filter with a constrained pole radius to decrease the roundoff-noise.

My thinking goes like this. We have some specification:

Maximum roundoff noise [dBFS]
Stopband attenuation. [dBFS]
Maximum Passband ripple [dB]
Passband upper frequency [fp |0 1| (Normalized)]
Stopband starting frequency [fs |fp 1| (Normalized)]

vs

computational cost.

Lets for an example use a 8:th order elliptic filter - but during testing we find that the variance in the roundoff noise is to high. One solution would be to change to higher prec. but the computational cost gets much higher - another solution might be to instead use a 12:th order filter with constrained pole radius that fits the specification - and such a solution will cost less than 50% extra computational power.

Graphical plots of a real result is upcoming - to be updated here.
Probably not the most confused programmer anymore on the XCORE forum.
Post Reply