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.
