|
|
| Line 1: |
Line 1: |
| = Introduction =
| |
|
| |
|
| Noise shaping is a technique that is often used in combination with
| |
| [[Dithering|dithering]] in the context of [[Quantize|quantization]]. The purpose of noise shaping
| |
| is to alter the spectral ''shape'' of the error that is introduced due
| |
| to [[Dithering|dithering]] and [[Quantize|quantization]]. Sometimes it may be desirable to
| |
| concentrate most of the noise power in a certain frequency band so
| |
| that the overall noise becomes less perceivable.
| |
|
| |
| = Applications =
| |
|
| |
| In the field of digital image processing a popular noise shaping
| |
| technique is known as "Floyd Steinberg dithering". This algorithm is
| |
| used to reduce the noise power at lower spatial frequencies which
| |
| results in an improved perceptual quality.
| |
|
| |
| For word length reductions of audio signals (i.e conversion from a 24bit
| |
| signal to a 16bit signal) noise shaping can be used to push quantization
| |
| noise out of perceptually important frequency bands to perceptually less
| |
| important frequency bands (>20 kHz) in order to reduce the perceived
| |
| loudness of the quantization noise.
| |
|
| |
| Another example is the
| |
| [[Super Audio Compact Disc|Super Audio CD]]. The SACD format represents
| |
| the most extreme trade-off between word length (1bit/sample) and
| |
| sampling frequency (2.8 Mhz). Due to the low word length high power
| |
| noise is introduced. With the help of noise shaping this noise is
| |
| usually rejected by about 120 dB in the audible band (0-20 kHz) at the
| |
| cost of a higher noise power in the ultrasonic range (+3.5dB for
| |
| frequencies above 100 kHz).
| |
|
| |
| = Limits & filter design restrictions =
| |
|
| |
| Noise Shaping without [[Dithering|''dithered'']] quantization is also possible but
| |
| less useful because the spectral characteristics of ''undithered''
| |
| quantization noise is not easily predictable and depends heavily on the
| |
| digital signal that is to be quantized.
| |
|
| |
| It's not possible to lower the noise power for ''all'' frequencies
| |
| via noise shaping because this would be in violation of the
| |
| rate/distortion principles. But it ''is'' possible to lower the noise
| |
| in certain frequency bands at the cost of increased noise power in other
| |
| frequency bands. This property is also known as the
| |
| ''noise shaping theorem''. It states that the average logarithmic magnitude response of
| |
| the corresponding noise shaping ''filter'' over the whole linear
| |
| frequency axis cannot drop below 0dB.
| |
|
| |
| The noise shaping application poses some restrictions on the fiter
| |
| design. The first sample of the filter's impulse response always equals
| |
| one. Also, interesting filters for noise shaping are so-called
| |
| ''minimum phase filters''. This type of filter performs optimal in the sense of
| |
| the noise shaping theorem because the average logarithmig magnitude
| |
| response ''is'' 0dB.
| |
|
| |
| For the special [[Super Audio Compact Disc|SACD]] case an additional
| |
| constraint comes into play: The magnitude response should not exceed
| |
| 3.5dB. This is enforced to prevent so-called limit cycles.
| |
|
| |
| The design of these kinds of filters is however not trivial. This is not
| |
| covered by this article.
| |
|
| |
| = Digital filter reminder =
| |
|
| |
| A common representation of digital filters is their ''transfer function''
| |
| in form of a rational polynomial in z^-1. These polynomials are defined
| |
| by their coefficients b0,b1,b2,...,a1,a2,...:
| |
|
| |
| <math>N(z) = \frac{B(z)}{A(z)} = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + \ldots + b_m z^{-m}}
| |
| {1 + a_1 z^{-1} + a_2 z^{-2} + \ldots + a_m z^{-m}} </math>
| |
|
| |
| where 'N' is the transfer function and 'm' the filter's ''order''.
| |
|
| |
| == FIR filters ==
| |
|
| |
| In case a1=a2=a3=...=0 the filter's impulse response is finite which is
| |
| why these filters are called FIR ([[Finite Impulse Response Filter|finite impulse response]])
| |
| filters as opposed to IIR
| |
| [[Infinite Impulse Response Filter|infinite impulse response]] filters. Note that in case of
| |
| FIR filters the coefficients b0, b1, ... correspond to the impulse
| |
| response. The impulse response is another way a filter can be
| |
| characterized and describes what the filter does to a unit pulse. Let's
| |
| assume for the moment that a1=a2=a3=...=0. The the impulse response
| |
| would be:
| |
|
| |
| x: 0 0 1 0 0 0 0 ... 0 0 0 ... (input signal)
| |
| y: 0 0 b0 b1 b2 b3 b4 ... 0 0 0 ... (output signal)
| |
|
| |
| It helps to think of negative powers of z in transfer functions as a
| |
| delay: z^-1 = delay of one sample. z^-2 = delay of two samples etc.
| |
| Since a digital filter is a linear and time-invariant system (LTI) it is
| |
| possible to determine the filter's output by a superposition of scaled
| |
| and shifted impulse responses. This leads to the following relationship
| |
| between input and output in the case of FIR filters:
| |
|
| |
| time domain Z-transform
| |
| ----------- -----------
| |
| y[k] = b0 x[k] + b1 x[k-1] + b2 x[k-2] + ... Y(z) = X(z) * B(z)
| |
|
| |
| where x[] is the input signal, y[] the output signal and k the current
| |
| sample index. This is a simple convolution directly computed using the
| |
| given impulse response. The convolution in the time domain corresponds
| |
| to a multiplication in the Z-domain.
| |
|
| |
| == IIR filters ==
| |
|
| |
| IIR filters are now a simple extension to the FIR case:
| |
|
| |
| time domain Z-transform
| |
| ----------- -----------
| |
| y[k] + a1 y[k-1] + a2 y[k-2] + ... Y(z) * A(Z) = X(Z) * B(Z)
| |
| = b0 x[k] + b1 x[k-1] + b2 x[k-2] + ...
| |
|
| |
| Solving for Y we get
| |
|
| |
| time domain Z-transform
| |
| ----------- -----------
| |
| y[k] = b0 x[k] + b1 x[k-1] + b2 x[k-2] + ... Y(Z) = X(Z) * N(z)
| |
| - a1 y[k-1] - a2 y[k-2] + ...
| |
|
| |
| which is the well known ''difference formula''. Its signal flow graph
| |
| representation is known as
| |
| [http://ccrma.stanford.edu/~jos/filters/Direct_Form_I.html direct form I].
| |
| Now, let's see how we can apply this kind of filter to quantization
| |
| errors...
| |
|
| |
| = Noise Shaper Implementation =
| |
|
| |
| Usually a filter's transfer function is called 'H'. In this article 'N'
| |
| has been used because this filter is also known as noise transfer
| |
| function (NTF) in the context of noise shaping. Later, another filter
| |
| 'H' is introduced that has a special relation to 'N'.
| |
|
| |
| == Direct Form I ==
| |
|
| |
| For convenience let us substitute the actual difference in the
| |
| difference formula for 'u'. Also, b0 must be equal to one for noise
| |
| shaping filters as mentioned above. This leads to
| |
|
| |
| y[k] = x[k] + u[k] with u[k] = b1 x[k-1] + b2 x[k-2] + ...
| |
| - a1 y[k-1] - a2 y[k-2] + ...
| |
|
| |
| Note that since the filter is supposed to be applied on the quantization
| |
| noise y[] corresponds to the filtered noise and x[] to the unfiltered
| |
| noise. Suppose s[] is the signal we're trying to quantize. The following
| |
| pseudo code snippet shows the inner loop over k of ''one possible''
| |
| noise shaping implementation:
| |
|
| |
| (1) compute u[k]
| |
| (2) wanted[k] = s[k] + u[k]
| |
| (3) quantized[k] = quantize( wanted[k] + dither() )
| |
| (4) y[k] = quantized[k] - s[k]
| |
| (5) x[k] = y[k] - u[k]
| |
|
| |
| That's about it. It may not be obvious at first, though. Take a closer
| |
| look at what happens at step 3. Note that the following is true
| |
|
| |
| time domain Z-transform
| |
| ----------- -----------
| |
| quantized[k] = s[k] + u[k] + x[k] Q(z) = S(z) + X(z) * N(z)
| |
| = s[k] + y[k] = S(z) + Y(z)
| |
|
| |
| u[k], wanted[k] are temporary variables. The values are not needed
| |
| outside this iteration, so you don't need to remember these. But we
| |
| ''do'' need to remember the past filtered and unfiltered error samples
| |
| y[] & x[] upto some amount of samples determined by the filter order.
| |
|
| |
| == Decoupling filtering and quantization ==
| |
|
| |
| [[Image:Noiseshaper3.png|thumb|right|200px]]
| |
| As mentioned, the previous example is ''one possible'' implementation.
| |
| It's not the smartest one. But it's hopefully easy to understand.
| |
| There's another common formulation. We now require some more transfer
| |
| function fiddeling. First, we add one and subtract one again:
| |
|
| |
| <math> N(z) = \frac{B(z)}{A(z)} = 1 - \frac{A(z)}{A(z)} + \frac{B(z)}{A(z)}</math>
| |
|
| |
| Since now the two fractions have a common denominator we can combine
| |
| them easily to one fraction. Also, let's factor z^-1 out of the
| |
| fraction:
| |
|
| |
| <math>N(z) = 1 - \frac{A(z)-B(z)}{A(z)} = 1 - z^{-1} \frac{z ( A(z)-B(z) )}{A(z)} =: 1 - z^{-1} H(z)</math>
| |
|
| |
| The last thing we did is writing H(z) for z*(A(Z)-B(Z)). Note that with
| |
| b0=1 it follows that z*(1-b0)=0. So, H(z) can be written as
| |
|
| |
| <math>H(z) = \frac{c0 + c1 z^{-1} + c2 z^{-2} + \ldots}{1 + a1 z^{-1} + a2 z^{-2} + \ldots} </math>
| |
|
| |
| with <math>c_{k-1} = a_k - b_k</math> for <math>k = 1 \ldots m</math>. H(z)
| |
| is some ''arbitrary'' filter. The big advantage now is that it is
| |
| possible to plug nearly any digital filter implementation for H(z) into
| |
| the noise shaping loop. One obvious approach is the difference formula
| |
| from above but there are other alternatives with desirable properties.
| |
| However, the purpose of this article is not to discuss digital filters
| |
| in general but their application as noise shaping filters. So,
| |
| throughout the rest of this article H(z) is simply treated as a black
| |
| box. This black box has a certain state and gives us an output sample
| |
| for each input sample. The modified pseudo code loop looks like this:
| |
|
| |
| (1) h_out[k] = H.filter(last_error)
| |
| (2) wanted[k] = s[k] - h_out[k]
| |
| (3) quantized[k] = quantize( wanted[k] + dither() )
| |
| (4) last_error = quantized[k] - wanted[k]
| |
|
| |
| which is equivalent to the shown signal flow graph where the unit
| |
| sample delay z^-1 is implemented by the state variable 'last_error':
| |
| Remember that this black box filter H(z) doesn't correspond the the
| |
| actual noise shaping filter N(z). Depending on the context published
| |
| coefficients for noise shaping correspond either to N(z), H(z), or
| |
| -H(z).
| |
|
| |
| = Example filters =
| |
|
| |
| The following two filters have been designed by Sebastian Gesemann
| |
| via non-linear optimization to match a slightly modified version of
| |
| the LAME MP3 encoder's ATH (absolute threshold of hearing) curve. They
| |
| were first published in the
| |
| [http://www.hydrogenaudio.org/forums/index.php?showtopic=47980&pid=551558&st=0 Hydrogenaudio forum]
| |
| including the noise transfer function's second order factorizations.
| |
|
| |
| == Sampling rate of 44100 Hz, derived from LAME's ATH formula ==
| |
|
| |
| The noise transfer function is given by
| |
| <math>N(z) = \frac{1 -1.1474 z^{-1} +0.5383 z^{-2} -0.3520 z^{-3} +0.3475 z^{-4}}
| |
| {1 +1.0587 z^{-1} +0.0676 z^{-2} -0.6054 z^{-3} -0.2738 z^{-4}}</math>
| |
|
| |
| The corresponding filter H is given by
| |
| <math>H(z) = \frac{2.2061 -0.4706 z^{-1} -0.2534 z^{-2} -0.6214 z^{-3}}
| |
| {1 +1.0587 z^{-1} +0.0676 z^{-2} -0.6054 z^{-3} -0.2738 z^{-4}}</math>
| |
|
| |
| == Sampling rate of 48000 Hz, derived from LAME's ATH formula ==
| |
|
| |
| The noise transfer function is given by
| |
| <math>N(z) = \frac{1 -1.3344 z^{-1} +0.7455 z^{-2} -0.4602 z^{-3} +0.3463 z^{-4}}
| |
| {1 +0.9030 z^{-1} +0.0116 z^{-2} -0.5853 z^{-3} -0.2571 z^{-4}}</math>
| |
|
| |
| The corresponding filter H is given by
| |
| <math>H(z) = \frac{2.2374 -0.7339 z^{-1} -0.1251 z^{-2} -0.6033 z^{-3}}
| |
| {1 +0.9030 z^{-1} +0.0116 z^{-2} -0.5853 z^{-3} -0.2571 z^{-4}}</math>
| |
