"High-End Audio" IIR filters

XCore Project reviews, ideas, videos and proposals.
User avatar
lilltroll
XCore Expert
Posts: 956
Joined: Fri Dec 11, 2009 3:53 am
Location: Sweden, Eskilstuna

"High-End Audio" IIR filters

Post by lilltroll »

Are you building an "HIFI" audio device including IIR filters that runs at fast samplingsfreq !?

Then you probably have problems with limit cycles, (even if you didn't know it)
http://en.wikipedia.org/wiki/Limit_cycle

Nonlinear !?, Well the rounding done for each sample introduces a non-linearity which is feedback into future samples, creating chaotic behaviour at least for a set of different filter coeffs. That makes the filter "sing" self sustained oscillations with a limited amplitude in the backgorund.
Typically bass management @ fs >=96 kHz creates a lot of chaotic energy decreasing the DNR Dynamic Noise Ratio of the filter.

The way to solve this is to use higher precision (high enough to make the error less than the quantization noise) in the calculation.
For an example using double float on a x86, uses 53 bits of digits for the internal states of the filter, and the multiply and add will be calculated with 80 bits (total) before rounding (http://en.wikipedia.org/wiki/X87) which is most often "good enough"

This filter-implementation uses 96 bits in the accumulator and 64 bits for the internal states.
What I have tested so far for realistic audio filters, the error out from the filter has equal magnitude as the quantization noise.

Compared to a software implementation with double, this implementation runs much faster.
You can run samples > 4.6 MHz/filter_order using one thread on a C5 device.

Check out:
https://www.xcore.com/projects/high-end ... ir-filters

PS. Using floating point doesn't make the problem with limit cycles to disappear, you must have enough bits in the significand (also coefficient or mantissa) DS.


User avatar
Gaston
Junior Member
Posts: 4
Joined: Wed Jan 27, 2010 10:39 am

Post by Gaston »

Very nice and solid work Lilltroll!

I guess there are multiple precision levels that can be used and they yield different results for different filters.
Is there any way to visualize the effect of increased precision and effect on DNR+THD maybe using noise and chirp as inputs?
Lets pick some easy and standard... Lets do a subwoofer crossover!
96kHz
First a LR4 (24 dB/octave Linkwitz-Riley) at 80Hz and maybe (I guess this could be implemented in the power amplifier?) a highpass BW3 (18 dB/octave ) at 14Hz for subsonic removal.

So what is needed to visualize, a couple of graphs for different inputs and for each graph a few traces of different color representing the precision?
User avatar
lilltroll
XCore Expert
Posts: 956
Joined: Fri Dec 11, 2009 3:53 am
Location: Sweden, Eskilstuna

Post by lilltroll »

It is a very good ideá to have some plots over the limit cycles, since there is no analytic solution to the problem, and I guess that the problem grows from almost nothing rather rapid over some threashold. A digital parametric EQ might work over a very wide freq.range (and Q-value as well), making it problem-free for higher freq. but starting to do limit cycles in the bass region.

Lets check out the freq. resp. of your example.
Attachments
IMG_08052012_181254.png
IMG_08052012_181254.png (15.14 KiB) Viewed 6868 times
IMG_08052012_181254.png
IMG_08052012_181254.png (15.14 KiB) Viewed 6868 times
User avatar
lilltroll
XCore Expert
Posts: 956
Joined: Fri Dec 11, 2009 3:53 am
Location: Sweden, Eskilstuna

Post by lilltroll »

First the filter coeff's in Q30.2 format

Code: Select all

B0={1073741824,-1072758407,0}
A0={1073741824,-1073741824,0}
B1={1072758407,-2145516815,1072758407}
A1={1073741824,-2146499330,1072758407}
B2={7332,14664,7332}
A2={1073741824,-2139532835,1065820340}
B3={7332,14664,7332}
A3={1073741824,-2139532835,1065820340}
First I need to calculate a very good reference.
That is done with help of vpa in MATLAB, using 128 decimals in base 10 (Direct Form I)

Code: Select all

for k=1:4096
     x=data(k,2);
     for i=1:4
         y =  vpa( (bd(i,1)*x + zp(i,1))/ad(i,1),DIGITS);
         zp(i,1) =  vpa(bd(i,2)*x + zp(i,2)- ad(i,2) * y,DIGITS);
         zp(i,2) =  vpa(bd(i,3)*x- ad(i,3) * y,DIGITS);
     x=y;
     end
     yp(k)=double(x);
 end
The input signal is a linear chirp sweept from DC to fs/2 over 4096 samples with the amplitude of 1.5E9
VPA128.png
(7.86 KiB) Not downloaded yet
VPA128.png
(7.86 KiB) Not downloaded yet
User avatar
Gaston
Junior Member
Posts: 4
Joined: Wed Jan 27, 2010 10:39 am

Post by Gaston »

That looks like the same frequency response as what I have:
Attachments
Subwoofer80HzLR4_14Hz_BW3.png
Subwoofer LP LR4 80Hz and HP BW3 14Hz
Subwoofer80HzLR4_14Hz_BW3.png (7.94 KiB) Viewed 6857 times
Subwoofer LP LR4 80Hz and HP BW3 14Hz
Subwoofer LP LR4 80Hz and HP BW3 14Hz
Subwoofer80HzLR4_14Hz_BW3.png (7.94 KiB) Viewed 6857 times
User avatar
lilltroll
XCore Expert
Posts: 956
Joined: Fri Dec 11, 2009 3:53 am
Location: Sweden, Eskilstuna

Post by lilltroll »

First, does my hocus pocus stuff work !?

Let's compare the reference result with the output from my XS1-L2.

Since the XMOS output is truncated to int32, the best possible result would be an an error of +-½

The result shows that the largest positive error is 0.536 and the largest negative error is -0.5278.
E.g. for some samples the result was not rounded to the correct value, BUT the error is only 1 LSB when it happens. Ex 0.536 should have been -0.464.
Attachments
VPA128vsXMOSext.png
VPA128vsXMOSext.png (7.06 KiB) Viewed 6846 times
VPA128vsXMOSext.png
VPA128vsXMOSext.png (7.06 KiB) Viewed 6846 times
User avatar
lilltroll
XCore Expert
Posts: 956
Joined: Fri Dec 11, 2009 3:53 am
Location: Sweden, Eskilstuna

Post by lilltroll »

Lets just test a DFII implementation with int64.
The problem with DFII is that the internal states need more than 32 bits!!

Code: Select all

MAX=0;
for k=1:4096
    x=int64(data(k,2));
    for i=1:4
        y =  (bd(i,1)*x + zp(i,1))/ad(i,1);
        zp(i,1) =  bd(i,2)*x + zp(i,2)- ad(i,2) * y;
        zp(i,2) =  bd(i,3)*x- ad(i,3) * y;
        MAX=max([MAX max(abs(zp))]);
    x=y;
    end
    accMAX(k)=MAX;
    MAX=0;
    ypi64(k)=double(x);
end
disp(MAX);
And the absolute error compared to the reference is not problem free.
Attachments
int64error.png
int64error.png (7.61 KiB) Viewed 6838 times
int64error.png
int64error.png (7.61 KiB) Viewed 6838 times
int64.png
int64.png (6.84 KiB) Viewed 6839 times
int64.png
int64.png (6.84 KiB) Viewed 6839 times
User avatar
lilltroll
XCore Expert
Posts: 956
Joined: Fri Dec 11, 2009 3:53 am
Location: Sweden, Eskilstuna

Post by lilltroll »

Let's change the implementation to Direct Form I, using a 64 bit accumulator and 32 bit internal states, (we now need 4pc. of 32 bits state, instead of 2 pc. 64 bits state )

Code: Select all

for k=1:4096
    x=int64(data(k,2));
    for i=1:4
        y = (bd(i,1)*x + bd(i,2)*x1(i) + bd(i,3)*x2(i)  - ad(i,2)*y1(i)-ad(i,3)*y2(i)) / ad(i,1);
        x2(i)=x1(i);
        x1(i)=x;
        y2(i)=y1(i);
        y1(i)=y;
        x=y;
    end
    accMAX(k)=max([max(abs([x2 x1 y1 y2]))]);
    MAX=0;
    ypd(k)=double(x);
end
Attachments
int64errorDFI.png
int64errorDFI.png (7.81 KiB) Viewed 6834 times
int64errorDFI.png
int64errorDFI.png (7.81 KiB) Viewed 6834 times
int64acc2.png
int64acc2.png (6.69 KiB) Viewed 6834 times
int64acc2.png
int64acc2.png (6.69 KiB) Viewed 6834 times
User avatar
lilltroll
XCore Expert
Posts: 956
Joined: Fri Dec 11, 2009 3:53 am
Location: Sweden, Eskilstuna

Post by lilltroll »

Lets take at look at single prec.

Code: Select all

for k=1:4096
    x=single(data(k,2));
    for i=1:4
        y = (bd(i,1)*x + bd(i,2)*x1(i) + bd(i,3)*x2(i)  - ad(i,2)*y1(i)-ad(i,3)*y2(i)) / ad(i,1);
        x2(i)=x1(i);
        x1(i)=x;
        y2(i)=y1(i);
        y1(i)=y;
        x=y;
    end
    accMAX(k)=max([max(abs([x2 x1 y1 y2]))]);
    MAX=0;
    yps(k)=double(x);
end
The error is so large that it can be seen direct in the output, partly due to the fact that the filter-coef's in truncated as well.
Attachments
singel error.png
singel error.png (7.12 KiB) Viewed 6832 times
singel error.png
singel error.png (7.12 KiB) Viewed 6832 times
User avatar
lilltroll
XCore Expert
Posts: 956
Joined: Fri Dec 11, 2009 3:53 am
Location: Sweden, Eskilstuna

Post by lilltroll »

Let's test x87 Extended double float (80-bits), with 64 bits double internal states.

Well we have a DNR of > 225 dB, but we would need a longer time-window to be sure of that.
Attachments
x87 error.png
x87 error.png (7.03 KiB) Viewed 6827 times
x87 error.png
x87 error.png (7.03 KiB) Viewed 6827 times
Post Reply