How to generate noise in frequency range with numpy?

Question

I have a main signal, for example sinus with period of 200 samples.

I would like to add a noise to this signal. The periods of "noise signal parts" should be in range for example 5-30 samples.

I thought that will be enough to generate multiple sinuses in this range with different randomly chosen amplitudes:

noise = np.sin(np.array(range(N))/0.7)*np.random.random(1) + np.sin(np.array(range(N))/1.1)*np.random.random(1) + np.sin(np.array(range(N))/1.5)*np.random.random(1) 

But this solution is still too much "deterministic" for my purpose.

How could I generate noise with randomly changing amplitude and period?

Answer 1

In MathWorks' File Exchange: fftnoise - generate noise with a specified power spectrum you find matlab code from Aslak Grinsted, creating noise with a specified power spectrum. It can easily be ported to python:

def fftnoise(f):
    f = np.array(f, dtype='complex')
    Np = (len(f) - 1) // 2
    phases = np.random.rand(Np) * 2 * np.pi
    phases = np.cos(phases) + 1j * np.sin(phases)
    f[1:Np+1] *= phases
    f[-1:-1-Np:-1] = np.conj(f[1:Np+1])
    return np.fft.ifft(f).real

You can use it for your case like this:

def band_limited_noise(min_freq, max_freq, samples=1024, samplerate=1):
    freqs = np.abs(np.fft.fftfreq(samples, 1/samplerate))
    f = np.zeros(samples)
    idx = np.where(np.logical_and(freqs>=min_freq, freqs<=max_freq))[0]
    f[idx] = 1
    return fftnoise(f)

Seems to work as far as I see. For listening to your freshly created noise:

from scipy.io import wavfile

x = band_limited_noise(200, 2000, 44100, 44100)
x = np.int16(x * (2**15 - 1))
wavfile.write("test.wav", 44100, x)

Answer 2

Instead of using multiple sinuses with different amplitudes, you should use them with random phases:

import numpy as np
from functools import reduce

def band_limited_noise(min_freq, max_freq, samples=44100, samplerate=44100):
    t = np.linspace(0, samples/samplerate, samples)
    freqs = np.arange(min_freq, max_freq+1, samples/samplerate)
    phases = np.random.rand(len(freqs))*2*np.pi
    signals = [np.sin(2*np.pi*freq*t + phase) for freq,phase in zip(freqs,phases)]
    signal = reduce(lambda a,b: a+b,signals)
    signal /= np.max(signal)
    return signal

Background: White noise means that the power spectrum contains every frequency, so if you want band limited noise you can add together every frequency within the band. The noisy part comes from the random phase. Because DFT is discrete, you only need to consider the discrete frequencies that actually occur given a sampling rate.