What's the relationship between the parameters of this noise generator and the gaussian curve it produces in the frequency domain?

Question

This algorithm produces a sequence $y$ having a magnitude spectrum that resembles a normal distribution, peaking at frequency $\omega$:

$$b[n]=(1-\sigma)b[n-1]+\frac{\sigma}{\omega} C$$ $$p[n]=p[n-1]+\omega(1+b[n])$$ $$y[n]=cos(p)$$

where $b[n]$ is a leaky integrator having feedback $1-\sigma$ and taking white noise $C$ as input, $p[n]$ is a phase accumulator having frequency $\omega$ modulated by brownian noise $b[n]$, and y[n] is the final signal, cosine at phase $p[n]$.

Python function:

def signal(f, w, n):
    """
    Arguments
    f - peak frequency in radians/sample
    w - bandwidth of frequency content
    n - length
    Returns
    sig - list of samples

    """
    br = 0
    p = 0
    sig = [0] * n
    for i in range(n):
        # br - leaky integrator having rate proportional to bandwidth,
        # integrating white noise at magnitude inversely proportional
        # to frequency.

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

        # p - phase increments at frequency f modulated by brownian noise.
        p = p + (1 + br) * f

        sig[i] = math.cos(p)

    return sig

To illustrate:

averaged magnitude spectrum at swept frequency from [0,π)

averaged magnitude spectrum with swept bandwidth from [0,0.5)

fourier transform of a sample sequence, hue indicates phase

wav

audio sample at 1760hz with $\sigma=0.005$

wav

averaged eight times

There are undesirable characteristics:

• The amplitude of the peak frequency has an unknown relationship to $\sigma$, the second animation shows the sharp curved response
• $\omega$ is not a perfect predictor of peak frequency, $\sigma$ allows noise to perturb the phase towards a range of $[0\omega,2\omega]$ from FM

Short of using optimisation and curve-fitting, how can the frequency content be predicted (e.g. by the gaussian function) and controlled?

Any additional insight appreciated

Answer 1

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.