import numpy as np

def noise(curve):
    # generate random phases for the signal
    phases = np.exp(1j * np.random.uniform(0.0, 2*np.pi, curve.shape))
    spectrum = curve * phases
    return np.fft.irfft(spectrum.real, norm='ortho')


def colored_noise(duration, rate, power=0.0):
    samples = int(duration * rate)
    freqs = np.fft.rfftfreq(samples, 1/rate)
    curve = np.zeros(samples // 2 + 1, dtype='float32')
    audible = (freqs >= 20) & (freqs < 20000)  # restrict to human audible range
    curve[audible] = freqs[audible] ** power  # compute curve
    curve *= np.sqrt(0.5 * rate / (curve**2).sum())  # normalize
    return noise(curve)


from scipy.io import wavfile

rate = 44100  # standard CD-quality sample rate
samples = colored_noise(10, rate, -0.5)
wavfile.write('noise.wav', rate, samples)


import sounddevice as sd

volume = 0.1

def callback(outdata, frames, time, status):
    outdata[:] = volume * colored_noise(frames / rate, rate, -1.0).reshape(-1, 1)


device = sd.query_devices(kind='output')
device_id = device['index']
rate = device['default_samplerate']
with sd.OutputStream(device=device_id, channels=1, blocksize=44100, callback=callback, samplerate=rate):
    print("Press enter to stop")
    input()
