How to implement filter poles in noise shaping for ΔΣ modulator? I am currently stuck trying to implement filter poles to stabilize a ΔΣ modulator. What I think I know is:

a noise shaping filter with only complex zeros along the unit circle destabilizes a 1 bit ΔΣ modulator when the order is beyond 2.
filter poles can be used to reduce the out-of-band noise, making the modulator more stable

For testing this approach, I am trying to realize a simple 2nd order modulator with 2 zeros and 2 poles.
Initially, I tried it without the poles and only with the two zeros which works and I implement it in the following way for testing (using python/numpy):
import numpy as np
roots = np.exp(1j*np.array([0.05, -0.05]))                  # two exemplary complex roots
b = np.poly(roots)                                          # and corresponding filter coeffs
errors = np.zeros_like(b)
signal = np.random.rand(100) - 0.5                          # any signal smaller than [-1,1]
samps = np.zeros_like(signal)
for iSAM in range(len(samps)):
    desired = signal[iSAM] + np.dot(b[1:],errors2[:-1])     # add the filtered errors to the input
    samps[iSAM] = np.sign(desired)                          # discretize to 1-bit
    errors = np.roll(errors, 1)
    errors[0] = samps[iSAM] - desired                       # calculate new error

As expected when FFT'ing the results, I see a frequency response like that of the desired filter.
Also as expected when adding another root at exp(i*0)=1 (creating a 3rd order modulator), the whole thing becomes unstable, even for small inputs.
I have no issue making filters with poles, and defining their frequency response. I am also able to take a time domain signal and filter it, with the designed biquad filter (2 zeros, 2 poles) and I obtain correct results.
The question is: how do I use the poles in the noise shaping algorithm. How would I even use a simple biquad filter for noise shaping ? I have watched many videos and papers today, but still don't understand it. Everyone draws these diagrams like this here from a video series i watched: https://youtu.be/IE8tU_10Hpg?t=63
A more detailed flow chart is in Stanley Lipshitz paper from 1991.
However, even if i follow this flow chart literally, I get no results. The state variables just explode. The guy in the video series also implemented a 3rd order filter with 3 zeros and 3 poles to define a stable 3rd order modulator. While I can reproduce the filter coefficients, and the Noise transfer function (NTF), I don't understand how to actually use the filter. In one video he makes a quick remark about testing the filter performance in a simulation. He says that performing the test is &quot;trivial&quot; or something like this. So obviously I am missing something very fundamental here :)
UPDATE:
After more days of trying, I am really at a loss. Can anyone please post an example of noise shaping with even a very simple IIR filter, even with only 2 zeros and 2 poles. The shaped noise spectrum should reproduce the magnitude response of the designed filter.
I am putting some bounty on it as much as I can. Because all the guides on the internet haven't helped me. I don't care about which programming language you use for it but the shaped spectrum should demonstrably produce the correct magnitude.
Best regards
I believe I have actually found an answer myself after even more research ;)
The crucial piece of information is the following: The filter function H(z) shown in many models is not the actual Noise Filter Function. Usually the noise filter should be some kind of high pass and I used this directly as H(z) erroneously.
Instead, the noise is shaped by the noise transfer function NTF(z) = 1-H(z).
Therefore, after designing the highpass noise transfer filter NTF(z), one needs to calculate the filter model via H(z) = 1-NTF(z).
It turns out that, if the NTF is a &quot;poleless&quot; filter (FIR filter) of form NTF(z) = 1 + Sum(b_n z^-n), then the filter function will be very similar: H(z) = 0 - Sum(b_n z^-n) where n starts from 1. The first coefficient will be 0.
I actually implemented this noise shaping directly in the past, by ignoring the leading coefficient 1 and inverting the sign of the filter result at the summing node with the signal. I came to this conclusion, because I picked it up somewhere but didn't understand its origins.
When the NTF contains poles, however, the H(z) will deviate more significantly. The denominator of NTF(z) and H(z) will be in fact identical. But the numerator will be X(z) = Sum(a_n z^-n) - Sum(b_n z^-n), where a_n are the pole coefficients, and b_n are the zero coefficients. As typically, the leading coefficients of the NTF(z) will be both 1 ( it is a minimum phase filter after all), the leading coefficient of H(z) will become 0. This produces the missing first coefficient naturally.
With this properly defined H(z), following the flow chart as given e.g. in Lipshitz 1991 paper, gives correct results. I used +/- 0.5 of the quantizer range as uniform white dither.
Update:
after correcting some minor mistakes, it all works as it should. In Python I use the Inverse Chebyshev Filter model to make the Noise Transfer Function. Then I correct the coefficients, to turn them into the H(z) filter model:
fDAC=20e6
b, a = signal.cheby2(5, 160, 20000, 'hp', fs=fDAC)
b /= b[0] # normalize according to Gerzon/Craven
b = a - b # turn NTF(z) into H(z) for noise shaping

Then the actual noise shaping (left some redundancy in the code, to keep it closer to Lipshitz 1991):
l = 200000
ts = np.arange(l)/fDAC
wav = 0.55 * np.sin(2*np.pi*3200*ts)
errbuf = np.zeros_like(b)
samps = np.zeros_like(wav)
for iSAM in range(l):
    desired = wav[iSAM] - np.dot(b, errbuf)
    samps[iSAM] = np.sign( desired )
    error = samps[iSAM] - desired 
    errbuf[0] = - error
    errbuf[0] = - np.dot(a, errbuf)
    errbuf = np.roll(errbuf, 1)

Note that no additive dither is being used. Even though the quantizer noise is not white itself, the noise shaping is efficient enough to remove the quantizer distortion. Undithered Delta Sigma modulators do have drawbacks tough, e.g. on stationary input, as can be read in numerous sources. If you want dither, it should be added right before quantization, i.e. inside the sign() function

Figure: Example of a 5th order delta sigma modulator. Vertical scale unit is full scale. Horizontal is Frequency in Hz. The orange curve is a designed 5th order inverse Chebyshev Highpass filter as Noise transfer function. Blue is the result with a ~-5dBFS sine wave.
PS: Since I originally set a bounty on this question, which was keeping me up at night: Please feel free to correct me and claim the bounty with a better answer of your own.