What's the relationship between the parameters of this noise generator and the gaussian curve it produces in the frequency domain?
This seems to work fine if everything is in rad/s. I think that the "problem" is with what you have implemented.

Here is a back of the envelope implementation of what you present:

Fs = 1000; % Sampling frequency in Hz
T = 4; % Total duration of the output signal
f = 125; % Centre frequency in Hz
s = 10; % Bandwidth in Hz

% Pre-allocate the output
y = zeros(1,T*Fs);

% Simple "buffers" for the previous sample in each filter
pb=0;
pp=0;

% For the duration of the signal...
for n = 0:(T*Fs)
    % Integrate b
    b = (1-2.*pi.*s./Fs).*pb + ((2.*pi.*s./Fs)./(2.*pi.*f./Fs)).*unifrnd(-1,1);
    % Integrate p
    p = pp + 2.*pi.*f./Fs.*(1+b);
    % Produce the output
    y(n+1) = cos(p);
    % Store the current values to the "buffers" to produce the n-1 indexing
    pb=b;
    pp=p;   
end;

% Frequency vector for the FFT assuming nbins=length(y)
f_v = Fs.*(0:(length(y)-1))./length(y);
% Spectrum of the produced signal
semilogy(f_v, abs(fft(y)));grid on; xlabel('Frequency (Hz)');ylabel('Amplitude');


For f=125, s=0:



Clean as a whistle, right on the f

For f=125, s=15.625



For f=125, s=62.5



Notice that I provide the bandwidth at fractions of f. The bandwidth cannot be wider than f.

I think that this is closer to the discrete Gaussian as it is more akin to a power law for very narrow bandwidths.

The deviation is at:


br = (1 - w ) * br + w * random.uniform(-1, 1) / max([0.001, f])


Why .../max([0.001,f])? The uniform() already returns numbers between $[-1..1]$ why do we need the additional "normalisation" there?



Overall, this is like an oscillator whose phase becomes a random variable. No matter how you are going to achieve it, this oscillator's phase must accumulate with an average rate of $\Delta \phi = \frac{2 \pi f}{Fs}$. The deviation from this would be proportional to the width around $\Delta \phi$. When setting a bandwidth, you are trying to limit the "phase growth" between two limits (but again, on average). So, for the Gaussian, the $\sigma$ (the standard deviation) does not correspond to the point in the curve that the output is at half the maximum. This occurs at $\approx \pm 1.1775 \sigma$, which you can apply as a correction factor to your s to modulate where you want your amplitude to be positioned approximately. But that would still depend on a few factors (for example the length of the signal and the whole sequence of random values that determines constructive and destructive interference throughout the signal).

Since you are in Python, you can construct such a "probabilistic" oscillator by drawing samples directly from a Gaussian $\mathcal{N}(\mu_f, \sigma_f)$ where $\mu_f, \sigma_f$ have frequency units as above) to obtain a series of $\Delta \phi$s and then simply add them up as you along (think of it as cumsum). Of course, this is not going to work if the hardware this oscillator is to run on can only generate uniform random numbers.

Hope this helps.