Okay, I admit to being a bit baffled at the moment....
It's probably the z-domain thingy. It takes a while to wrap one's head around it.
I looked up noise shaping in Wikipedia and foundy(n) = x(n)+A.E(x(n-1))+B.E(x(n-2))+C.E(x(n-3))+D.E(x(n-4))+
E.E(x(n-5))+F.E(x(n-6))+G.E(x(n-7))+H.E(x(n-8))+I.E(x(n-9))
I also found some code which was using a a filter with 9 coefficients so I implemented the noise shaping in lossyWAV like that, i.e. output = input - coeff[0..8] x quantization_error[0..-8], where quantization_error = output - input.
The 1st problem with this wikipedia article is that it's not really obvious what E is. Is it the unfiltered or the filtered noise? Btw, output-input isn't the the quantization error. It's the already-filtered error. So, in your case you'll get an all-pole filter which is a totally different beast than a FIR filter where the actual quantization errors are used. The difference is subtle: Note, that in the picture I made I pick up the signal after the feedback and right before dither and quantization noise is added to compute the "quantization error" (unfiltered noise).
The 2nd problem with this wikipedia article is that it doesn't say anything about FIR or IIR filters and whether and/or how they can be used and what type of filter is actually described there.
That said: Regardless of what E is, their noise shaping formula is equivalent to the structure I drew where H(z) either corresponds to an all-pole-IIR or a FIR filter.
Initially, this kept crashing until I divided all coefficients by coeff[0] and then disregarded coeff[0] (per David's code).
By "crashing" I guess you mean the filter went unstable. You probably used the coefficients in a wrong way. It might be a sign problem (sign of E is wrong) or you got the wrong E (filtered noise instead of unfiltered noise).
Looking again at the Wikipedia article, it appears that I have omitted to include dither in the calculation.
I still feel as if I'm groping in the dark here and would gratefully accept any advice, pointers, etc.
If you omit dither you can't guarantee the quantization error to be white/uncorrelated. The noise shaping stuff still works but you may get unexpected results because the filter is supposed to be applied on white/uncorrelated noise. So, that's why at least rectangular dithering should be used.
More explanations and pseudo code following...
You might have missed some informations regarding the picture I drew
X : input signal
Y : output signal ( = input + filtered noise )
E : dither & quantization noise (unfiltered white noise please)
+ : is obviously mixing two signals. Note it can also be used
for subtraction (source line(s) marked with a minus)
Also, quantization is modelled as mixing the signal with
errors.
[z^-1] : This is a simple filter: A delay of one sample
[H(z)] : This is any filter you like to use
So, suppose you have some given filter coefficients for H(z):
b[] = {b[0],b[1],b[2],...,b[n]}; // array, indexed starting at 0
a[] = { 1 ,a[1],a[2],...,a[m]}; // array, we don't need a[0]
The index actually corresponds to the power of 1/z for the z-domain
interpretation, 'b' holds the numerator coefficients and 'a' holds
the denominator coefficients.
x[k] and y[k] are the input and output samples.
We also need some filter memory with exactly max(n+1,m) samples. Let's
write fifo[0] for the last sample we added to the fifo, fifo[1] was
the last sample in the previous loop and so on...
Then, the inner loop over 'k' would look like this:
{
wanted_temp = x[k] + fifo[0] * b[0]
+ fifo[1] * b[1]
+ ..............
+ fifo[n] * b[n];
y[k] = quantize( wanted_temp + dither );
qerror_temp = wanted_temp - y[k];
new_fifo_sample = qerror_temp - fifo[0] * a[1]
- fifo[1] * a[2]
- ..............
- fifo[m-1] * a[m];
fifo_add( fifo, new_fifo_sample );
// Now: fifo[0] == new_fifo_sample
}
For implementing H(z) I used the direct form II structure where the delay-line is shared among the recursive and non-recursive filter parts.
The 4th order filter I was suggesting for 24->16 bit word length reduction @ 44 kHz sampling frequency leads to the following coefficients for H(z):
b[0..3] = { 2.2061 , -0.4707 , -0.2534 , -0.6213 };
a[1..4] = { 1.0587 , 0.0676 , -0.6054 , -0.2738 };
Again: H(z) is NOT the transfer function of the noise shaper, it is G(z) = 1 - z^-1 * H(z).
Note: This post comes with no warrenty and might contain errors.
Cheers,
SG
Topic: lossyWAV Development (Read 573711 times)
previous topic - next topic