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):
    f - peak frequency in radians/sample
    w - bandwidth of frequency content
    n - length
    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


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


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

  • $\begingroup$ Ideally, how do you expect this to work? Say for instance, you expect to set something like $\sigma=100, \omega=2 \pi 440 rad/s$ and this produces $440 Hz \pm \approx 200$ at full width at half maximum? What I am trying to say is that it may simply be a matter of "fixing" your units and "catering for your gains" (for example, you have there an $h = \left[1-\sigma, \frac{\sigma}{\omega} \right]$ whose "gain" changes as a function of both $\sigma$ and $\omega$. Are you working off of a specific reference for this? $\endgroup$ – A_A Jun 26 '18 at 15:16
  • $\begingroup$ Thanks for reading and replying, seems like you've got the right idea, though it's not clear how to fix the units to create a specific width like you say (I wasn't aware of FWHM as a concept but it fits the problem perfectly). Eg at first the spectrum widened with increasing $\omega$ until the noise was made inversely proportional. The gain of the peak frequency varies with $\sigma$. There's no specific reference for the whole system no, just the definition of FM combined with a brownian walk as modulator. Maybe I should plot the current relationship between $\sigma$ to amplitude & FWHM? $\endgroup$ – JavariG Jun 27 '18 at 3:52

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

% 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

% 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:

enter image description here

Clean as a whistle, right on the f

For f=125, s=15.625

enter image description here

For f=125, s=62.5

enter image description here

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.