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.
"High-End Audio" IIR filters
-
- XCore Expert
- Posts: 956
- Joined: Fri Dec 11, 2009 3:53 am
- Location: Sweden, Eskilstuna
-
- Junior Member
- Posts: 4
- Joined: Wed Jan 27, 2010 10:39 am
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?
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?
-
- XCore Expert
- Posts: 956
- Joined: Fri Dec 11, 2009 3:53 am
- Location: Sweden, Eskilstuna
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.
Lets check out the freq. resp. of your example.
You do not have the required permissions to view the files attached to this post.
-
- XCore Expert
- Posts: 956
- Joined: Fri Dec 11, 2009 3:53 am
- Location: Sweden, Eskilstuna
First the filter coeff's in Q30.2 format
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)
The input signal is a linear chirp sweept from DC to fs/2 over 4096 samples with the amplitude of 1.5E9
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}
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
You do not have the required permissions to view the files attached to this post.
-
- Junior Member
- Posts: 4
- Joined: Wed Jan 27, 2010 10:39 am
That looks like the same frequency response as what I have:
You do not have the required permissions to view the files attached to this post.
-
- XCore Expert
- Posts: 956
- Joined: Fri Dec 11, 2009 3:53 am
- Location: Sweden, Eskilstuna
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.
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.
You do not have the required permissions to view the files attached to this post.
-
- XCore Expert
- Posts: 956
- Joined: Fri Dec 11, 2009 3:53 am
- Location: Sweden, Eskilstuna
Lets just test a DFII implementation with int64.
The problem with DFII is that the internal states need more than 32 bits!!
And the absolute error compared to the reference is not problem free.
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);
You do not have the required permissions to view the files attached to this post.
-
- XCore Expert
- Posts: 956
- Joined: Fri Dec 11, 2009 3:53 am
- Location: Sweden, Eskilstuna
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
You do not have the required permissions to view the files attached to this post.
-
- XCore Expert
- Posts: 956
- Joined: Fri Dec 11, 2009 3:53 am
- Location: Sweden, Eskilstuna
Lets take at look at single prec.
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.
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
You do not have the required permissions to view the files attached to this post.
-
- XCore Expert
- Posts: 956
- Joined: Fri Dec 11, 2009 3:53 am
- Location: Sweden, Eskilstuna
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.
Well we have a DNR of > 225 dB, but we would need a longer time-window to be sure of that.
You do not have the required permissions to view the files attached to this post.