# limit to possible noise shaping?

I want to do noise shaping in a 100kHz, 16 bit application, so as to shift all quantization noise to the 25khz-50kHz band, with minimal noise in the DC-25kHz band.

I set up mathematica to create a 31-sample error filter kernel via reinforcement learning which works well: After a little learning, I can get about ~16dB enhancement of high frequency noise for about the same amount of reduction in the low frequency band (the central line is the unshaped dither noise level). This is in line with the "Gerzon-Craven" noise shaping theorem.

Now to my problem:

I cannot manage to shape the noise even more even after extensive learning, although the Gerzon-Craven theorem doesnt forbid it. For example, it should be possible to achieve 40 dB reduction in the low band and 40 dB enhancement in the high band.

So is there another fundamental limit I am running into ?

I tried looking into the Shannon noise/sampling/information theorems but after fiddling with that a while, I only managed to derive a single limit from it: the Gerzon-Craven theorem, which seems to be a direct outcome of the Shannon theorem.

Any help is appreciated.

First off, the filter kernel that produces the above noise shaping, Note that the most recent sample is on the right side. Numeric values of the BarChart rounded to .01: {-0.16, 0.51, -0.74, 0.52, -0.04, -0.25, 0.22, -0.11, -0.02, 0.31, -0.56, 0.45, -0.13, 0.04, -0.14, 0.12, -0.06, 0.19, -0.22, -0.15, 0.4, 0.01, -0.41, -0.1, 0.84, -0.42, -0.81, 0.91, 0.75, -2.37, 2.29} (Not exactly the bar char but produces similar curve)

Another note about error feedback implementation:

I tried two different implementations of the error feedback. First I compared the rounded output sample to the desired value and used this deviation as error. Second I compared the rounded output sample to the (input+error feedback). Although both methods produce quite different kernels, both seem to level off at about the same noise shaping intensity. The data posted here uses the second implementation.

Here is the code that is used to calculate the digitized wave samples. step is the stepsize for rounding. wave is the undigitized waveform (typically just zeros when no signal is applied).

TestWave[kernel_?VectorQ] :=
Module[{k = kernel, nf, dith, signals, twave, deltas},
nf = Length@k;
dith = RandomVariate[TriangularDistribution[{-1, 1}*step], l];
signals = deltas = Table[0, {l}];
twave = wave;
Do[
twave[[i]] -= k.PadLeft[deltas[[;; i - 1]], nf];
signals[[i]] = Round[twave[[i]] + dith[[i]], step];
deltas[[i]] = signals[[i]] - twave[[i]];
, {i, l}];
signals
]


The reinforcement method:

The "score" is computed by looking at the noise power spectrum. The goal is to minimize noise power in the band DC-25kHz. I am not penalizing noise in the high frequency band, so arbitrarily high noise there would not diminish score. I am introducing noise in the kernel weights to learn. Maybe, therefore, I am in a (very broad and deep) local minimum, but I consider this extremely unlikely.

Comparison to standard filter design:

Mathematica allows to generate filters iteratively. These can have much better contrast than 36 dB when their frequency response is plotted; up to 80-100 dB. Numeric values: {0.024, -0.061, -0.048, 0.38, -0.36, -0.808, 2.09, -0.331, -4.796, 6.142, 3.918, -17.773, 11.245, 30.613, -87.072, 113.676, -87.072, 30.613, 11.245, -17.773, 3.918, 6.142, -4.796, -0.331, 2.09, -0.808, -0.36, 0.38, -0.048, -0.061, 0.024}

However, when applying those in the actual noise shaping, they are (a) clamped to the same ~40dB contrast, (b) perform worse than the learned filter doing actually no noise attenuation.

• +1, very interesting question. Have you tried increasing the filter's order above 31 taps? 40dB suppression sounds a bit high for a 31 tap FIR. – A_A Feb 15 '18 at 10:06
• @Olli, I don't believe I understand completely. I can post the filter kernel if that is what you are interested in. In blunt words, there are oscillatory weights which forces the error to alternative -> shifts it to high frequencies. – tobalt Feb 15 '18 at 10:08
• @tobalt from "classical" filter design, it's an expected result that longer filters are steeper and/or have more attenuation in the stop band and/or have less ripple in the pass band. Now, my guess is that your reinforcement method rewards steepness more than attenuation after some point; what's the method that you use to reinforce? – Marcus Müller Feb 15 '18 at 10:25
• You might want to have a look at the Filter Design section of Mathematica. Perhaps you can simply define the specifications of your filter and use one of the existing techniques to return a filter that satisfies them. – A_A Feb 15 '18 at 11:33
• Then that is definitely (optionally iterative) filter design. Get your filter specifications (exactly as you posted them here) and try to create a filter via this function (the simplest of its type) and see what it comes up with. It would be nice to see the coefficients that function comes up with versus the ones the re-inforcement learning comes back with. Also, note the sort of filter order it comes up with, I am guessing it would be higher than 31. Does it have to be "adaptive" to the signal by the way? – A_A Feb 15 '18 at 11:46

# Basic dithering without noise shaping

Basic dithered quantization without noise shaping works like this:

Figure 1. Basic dithered quantization system diagram. Noise is zero-mean triangular dither with a maximum absolute value of 1. Rounding is to nearest integer. Residual error is the difference between output and input, and is calculated for analysis only.

The triangular dither increases the variance of the resulting residual error by a factor of 3 (from $\frac{1}{12}$ to $\frac{1}{4}$) but decouples the mean and variance of the net quantization error from the value of the input signal. That means that the net error signal is uncorrelated with the input but higher moments are not decoupled, so it's not truly completely independent random error, but no one has determined that people can hear any dependency of the higher moments in the net error signal on the input signal in an audio application.

With independent additive residual error we would have a simpler model of the system:

Figure 2. Approximation of basic dithered quantization. Residual error is white noise.

In the approximate model the output is simply input plus independent white noise residual error.

# Dithering with noise shaping

I can't read Mathematica very well so instead of your system I'll analyze the system from Lipshitz et al. "Minimally audible noise shaping" J. Audio Eng. Soc., Vol.39, No.11, November 1991:

Figure 3. Lipshitz et al. 1991 system diagram (adapted from their Fig. 1). Filter (italicized in the text) includes in it a one sample delay so that it can be used as an error feedback filter. Noise is triangular dither.

If the residual error is independent from current and past values of signal A, we have a simpler system:

Figure 4. An approximate model of the Lipshitz et al. 1991 system. Filter is the same as in Fig. 3 and includes in it a one sample delay. It is no longer used as a feedback filter. Residual error is white noise.

In this answer I will work with the more easily analyzed approximate model (Fig. 4). In the original Lipshitz et al. 1991 system, Filter has a generic infinite impulse response (IIR) filter form that covers both IIR and finite impulse response (FIR) filters. In the following we will assume that Filter is a FIR filter, as I believe based on my experiments with your coefficients that that is what you have in your system. The transfer function of Filter is:

$$H_\mathit{Filter}(z) = -b_1z^{-1} - b_2z^{-2} - b_3z^{-3} - \ldots$$

The factor $z^{-1}$ represents a one-sample delay. In the approximate model there is also a direct summation path to output from residual error. This gets summed with the negated output of Filter, forming the full noise shaping filter transfer function:

$$H(z) = 1 - H_\mathit{Filter}(z) = 1 + b_1z^{-1} + b_2z^{-2} + b_3z^{-3} + \ldots.$$

To go from your Filter coefficients, which you list in order $\ldots, -b_3, -b_2, -b_1$, to the full noise shaping filter transfer function polynomial coefficients $1, b_1, b_2, b_3, \ldots$, the sign of the coefficients is changed to account for the negation of Filter output in the system diagram, and the coefficient $b_0 = 1$ is appended to the end (by horzcat in the Octave script below), and finally the list is reversed (by flip):

pkg load signal
b = [-0.16, 0.51, -0.74, 0.52, -0.04, -0.25, 0.22, -0.11, -0.02, 0.31, -0.56, 0.45, -0.13, 0.04, -0.14, 0.12, -0.06, 0.19, -0.22, -0.15, 0.4, 0.01, -0.41, -0.1, 0.84, -0.42, -0.81, 0.91, 0.75, -2.37, 2.29];
c = flip(horzcat(-b, 1));
freqz(c)
zplane(c)


The script plots the magnitude frequency response and the zero locations of the full noise shaping filter:

Figure 5. Magnitude frequency response of the full noise-shaping filter.

Figure 6. Z-plane plot of poles ($\times$) and zeros ($\circ$) of the filter. All the zeros are inside the unit circle, so the full noise-shaping filter is minimum-phase.

I think the problem of finding the filter coefficients can be reformulated as the problem of designing a minimum-phase filter with a leading coefficient of 1. If there are inherent limitations to the frequency response of such filters, then these limitations are transferred to equivalent limitations in noise shaping that uses such filters.

# Conversion from all-pole design to minimum-phase FIR

A procedure for design of different but in many ways equivalent filters are described in Stojanović et al., "All-Pole Recursive Digital Filters Design Based on Ultraspherical Polynomials", Radioengineering, vol 23, no 3, September 2014. They calculate denominator coefficients of the transfer function of an IIR all-pole low-pass filter. Those always have a leading denominator coefficient of 1 and have all poles inside the unit circle, a requirement of stable IIR filters. If those coefficients are used as the coefficients of the minimum-phase FIR noise shaping filter, they will give an inverted high-pass frequency response compared to the low-pass IIR filter (transfer function denominator coefficients becoming numerator coefficients). In your notation one set of coefficients from that article is {-0.0076120, 0.0960380, -0.5454670, 1.8298040, -3.9884220, 5.8308660, -5.6495140, 3.3816780}, which could be tested for the noise shaping application although it isn't exactly to the specification:

Figure 7. Magnitude frequency response of the FIR filter using coefficients from Stojanović et al. 2014.

Figure 8. Pole-zero plot of the FIR filter using coefficients from Stojanović et al. 2014.

The all-pole transfer function is:

$$H(z) = \frac{1}{1 + a_1z^{-1} + a_2z^{-2} + a_3z^{-3} + \ldots}$$

So, you can design a stable all-pole IIR low-pass filter and use the $a$ coefficients as $b$ coefficients to get a minimum-phase high-pass FIR filter with a leading coefficient of 1.

To design an all-pole filter and to convert that to a minimum-phase FIR filter, you will not be able to use IIR filter design methods that start from an analog prototype filter and map the poles and zeros into the digital domain using bilinear transform. That includes cheby1, cheby2, and ellip in Octave and Python's SciPy. These methods will give zeros away from z-plane origin so the filter will not be of the required all-pole type.

# Answer to the theoretical question

If you don't care how much noise there will be at frequencies above quarter of the sampling frequency, then Lipshitz et al. 1991 addresses your question directly:

For such weighting functions, which go to zero over part of the band, there is no theoretical limit to the weighted noise-power reduction obtainable from the circuit of Fig. 1. This would be the case if, for example, one assumes that the ear has zero sensitivity between, say, 20 kHz and the Nyquist Frequency, and chooses the weighting function to reflect this fact.

Their Fig 1. shows a noise shaper with a generic IIR filter structure with both poles and zeros, so different to the FIR structure that you have at the moment, but what they say applies also to that, because a FIR filter impulse response can be made arbitrarily close to the impulse response of any given stable IIR filter.

# Octave script for filter design

Here is an Octave script for coefficient calculation by another method that I think is equivalent to the Stojanovici et al. 2014 method parameterized as $\nu=0$ with the right choice of my dip parameter.

pkg load signal
N = 14; #number of taps including leading tap with coefficient 1
att = 97.5; #dB attenuation of Dolph-Chebyshev window, must be positive
dip = 2; #spectrum lift-up multiplier, must be above 1
c = chebwin(N, att);
c = conv(c, c);
c /= sum(c);
c(N) += dip*10^(-att/10);
r = roots(c);
j = (abs(r(:)) <= 1);
r = r(j);
c = real(poly(r));
c .*= (-1).^(0:(N-1)); #if this complains, then root finding has probably failed
freqz(c)
zplane(c)
printf('%f, ', flip(-c(2:end))), printf('\n'); #tobalt's format


It starts with a Dolph-Chebyshev window as the coefficients, convolves it with itself to double the transfer function zeros, adds to the middle tap a number that "lifts up" the frequency response (considering the middle tap as being at zero time) so that it is everywhere positive, finds the zeros, removes zeros that are outside the unit circle, converts the zeros back to coefficients (the leading coefficient from poly is always 1), and flips the sign of every second coefficient to make the filter high-pass. Results from (an older but almost equivalent version of) the script look promising:

Figure 9. Magnitude frequency response of the filter from (an older but almost equivalent version of) the above script.

Figure 10. Pole-zero plot of the filter from (an older but almost equivalent version of) the above script.

The coefficients from (an older but almost equivalent version of) the above script in your notation: {0.357662, -2.588396, 9.931419, -26.205448, 52.450624, -83.531276, 108.508775, -116.272581, 102.875781, -74.473956, 43.140431, -19.131434, 5.923468}. The numbers are large which could lead to numerical problems.

# Octave implementation of noise shaping

Finally, I did my own implementation of noise shaping in Octave and don't get problems like you did. Based on our discussion in comments, I think the limitation in your implementation was that the noise spectrum was evaluated using a rectangular window a.k.a. "no windowing", which spilled the high-frequency spectrum to the low frequencies.

pkg load signal
N = length(c);
M = 16384; #signal length
input = zeros(M, 1);#sin(0.01*(1:M))*127;
er = zeros(M, 1);
output = zeros(M, 1);
for i = 1:M
A = input(i) + er(i);
output(i) = round(A + rand() - rand());
for j = 2:N
if (i + j - 1 <= M)
er(i + j - 1) += (output(i) - A)*c(j);
endif
endfor
endfor
pwelch(output, max(nuttallwin(1024), 0), 'semilogy');


Figure 11. Quantization noise spectral analysis from the above Octave implementation of noise shaping for constant zero input signal. Horizontal axis: Normalized frequency. Black: no noise shaping (c = [1];), red: your original filter, blue: the filter from section "Octave script for filter design".

Figure 12. Time domain output from the above Octave implementation of noise shaping for constant zero input signal. Horizontal axis: sample number, vertical axis: sample value. Red: your original filter, blue: the filter from section "Octave script for filter design".

The more extreme noise shaping filter (blue) results in very large quantized output sample values even for zero input.

• @MattL. I thought wrongly at first that tobalt has an all-pole filter. I rewrote my answer when I realized it is a FIR filter with the first coefficient 1. Also Gerzon-Craven is reported to say that the filter needs to be minimum phase to be optimal, and tobalt's optimized coefficients also give a minimum phase filter. Those requirements are equivalent to what the coefficients of IIR all-pole filters have so I suggest borrowing design methods from there. A standard IIR would be an option too. – Olli Niemitalo Feb 16 '18 at 12:17
• I have isolated the error: My implementation produces the same waveform (in time) as yours does. However the Abs[Fourier[wave]] function seems to run into some internal overflow/underflow, because the returned spectrum looks different (higher floor) – tobalt Feb 18 '18 at 10:42
• @Olli Niemitalo Ok it seems like the FFT in octave uses automatic windowing possibly ? After applying a Hann window to the waveform I can get "correct" FFTs. I will briefly test the integrity of this approach and eventually continue learning and post the outcome. Thanks for all your efforts. I have marked your post as answer. – tobalt Feb 18 '18 at 11:15
• @robertbristow-johnson I think it's all consistent as it is. I removed an equation where H(z) was for a recursive filter with 1 as the numerator. But it's a FIR filter in tobalt's case. I suspect you may think that it becomes a recursive filter because there is a feedback loop. But dithered quantization is in the loop doing its thing cutting the path from the filter output to the residual. – Olli Niemitalo Feb 18 '18 at 18:27
• Also Lipshitz et al. 1991 use $a$ and $b$ with the opposite meanings, a practice I was schooled out of here on dsp.stackexchange.com for being non-standard. – Olli Niemitalo Feb 18 '18 at 18:37