Wilbert's website at SocSci

> Software> Python> Noise

python/noise.html 2019-10-30

Noise

White noise

Making noise in python is very simple. To make 1000 gaussian white noise samples do:

#!/usr/bin/env python3
import numpy as np
n = 1000
np.random.normal(0.0,  1.0, n)

These samples are independant and have a gaussian distribution with mean = 0.0 and SD = 1.0:

Note that for noise to be white there is absolutely no requirement to have a gaussian distribution. The following yields (non gaussian) white noise:

np.random.randint(0, 2, n)

Pink noise

Pink noise requires making a spectrum with power spectral density (PSD) S = 1/f and reverse fourier transforming it. Its quite simple in Python:

#!/usr/bin/python3
import numpy as np

def spectrum_noise(spectrum_func, samples=1024, rate=44100):
	""" 
	make noise with a certain spectral density
	"""
	freqs = np.fft.rfftfreq(samples, 1.0/rate)            # real-fft frequencies (not the negative ones)
	spectrum = np.zeros_like(freqs, dtype='complex')      # make complex numbers for spectrum
	spectrum[1:] = spectrum_func(freqs[1:])               # get spectrum amplitude for all frequencies except f=0
	phases = np.random.uniform(0, 2*np.pi, len(freqs)-1)  # random phases for all frequencies except f=0
	spectrum[1:] *= np.exp(1j*phases)                     # apply random phases
	return np.fft.irfft(spectrum)                         # return the reverse fourier transform

def pink_spectrum(f, f_min=0, f_max=np.inf, att=np.log10(2.0)*10):
	"""
	Define a pink (1/f) spectrum
		f     = array of frequencies
		f_min = minimum frequency for band pass
		f_max = maximum frequency for band pass
		att   = attenuation per factor two in frequency in decibel.
		        Default is such that a factor two in frequency increase gives a factor two in power attenuation.
	"""
	# numbers in the equation below explained:
	#  0.5: take the square root of the power spectrum so that we get an amplitude (field) spectrum 
	# 10.0: convert attenuation from decibel to bel
	#  2.0: frequency factor for which the attenuation is given (octave)
	s = f**-( 0.5 * (att/10.0) / np.log10(2.0) )  # apply attenuation
	s[np.logical_or(f < f_min, f > f_max)] = 0    # apply band pass
	return s

The code should be pretty self explenatory. The default attenuation gives standard 1/f pink noise. This is sometimes called 3 dB rolloff.

Usage

To make a wav-file with this noise in the 5 - 10 kHz band scaled to the maximum volume ofg wav-files, use:

from scipy.io import wavfile
x = spectrum_noise(lambda x:pink_spectrum(x, 5000, 10000), 4096)
x = np.int16(x/ max(abs(x)) * (2**15 - 1)) 
wavfile.write("pink.wav", 44100, x)

Modulation

For 44100 samples of high pass noise above 140 Hz with 6 dB rolloff and 2 Hz modulation with random modulation phase, use:

from scipy.io import wavfile
samples = 44100 # number of samples
rate = 44100    # rate (Hz)
f_mod = 2       # modulation frequency (Hz)
m = 1           # modulation depth
x = spectrum_noise(lambda x:pink_spectrum(x, 140, np.inf, 6.0), samples, rate)
x *= 1.0 + m * np.sin(2*np.pi * ( f_mod * (np.arange(0, samples)/rate + np.random.rand()) )) # modulation
factor = (2**13-1) / np.std(x)  # scale rms to 2^15 / 4 \approx -12 dBFS RMS
#factor = (2**15-1) / max(abs(x)) # scale peak to 2^15 = 0 dBFS peak
x = np.int16(x * factor)  # scale and convert to 16 bit integer


wavfile.write("pink.wav", rate, x)

Finding the level

We now have a file nicely normalised to -12dBFS RMS. But how does this compare to the sound that this noise has to be added to? Plenty of programs that can find an audio file's RMS:

normalize-audio -n sound.wav 
ffmpeg -i sound.wav -filter:a volumedetect -f null /dev/null
ffmpeg-normalize -n sound.wav -f --print-stats -nt rms

Note that there are other way of determining the level of a file: peak and ebu.

Normalisation (scaling)

ffmpeg-normalize -nt ebu -t -12  -ofmt wav -f pink.wav -o pink2.wav

Mixing

Now that we know how loud the noise has to become, we can add it to the sound:

sox -m pink.wav sound.wav pink3.wav
sox pink.wav pink2.wav --plot gnuplot -V2 loudness  | gnuplot -p

Misc commands

cvlc --play-and-exit pink.wav # play with cvlc
python3 -c "from scipy.io import wavfile; print(wavfile.read('pink.wav')[1].shape)" # show number of samples