7
$\begingroup$

How can I generate arbitrarily high quality $1/f$ noise, for use in a model?

My model involves a lot of feedback, over a large number of iterations, with a very high bandwidth, so I'd like the $1/f$ noise to be as "ideal" as possible.

Assuming there's not an ideal solution, it would be nice to have an approach that allowed "noise quality" to be a parameter, so I could generate arbitrarily good noise.

Using scipy/numpy right now for exploring the model, but I'll probably be reimplementing in C.

Improve this question
$\endgroup$
2
  • $\begingroup$ FYI the accepted answer here is incorrect. See my comment below. $\endgroup$ – DanielSank Jan 28 '16 at 2:22
  • $\begingroup$ Where is 'n' used in these set of formulas?<br> What is the unit of T? It should be Hz if f<sub>k</sub> is a frequency, not Time. <br> How to determine k<sub>max</sub>? I want to generate pink noise between lowest freuency and 44 kHz.<br> And pink noiise between 1 and 2 kHz.<br> Which values to use? $\endgroup$ – Peter Feb 26 '20 at 10:05
9
$\begingroup$

You can generate a noise sequence with whatever noise spectrum you want (including $1/f$, also known as pink noise) by generating the noise coefficients in spectral space. The magnitudes of the coefficients should be chosen to give the desired spectrum and the phases should be chosen randomly. You then simply perform an inverse Fourier transform to give the sequence values. The following algorithm will generate 1/f noise:

  1. Determine the number of points, $n$, and the length, $T$, of your sequence. This also determines the spectral space with wave numbers $-k_{\max}$ to $k_{\max}$ which correspond to frequency values $f_k=kT/(2\pi)$.
  2. Set the magnitudes of your spectral coefficients: $C_k = 1/\left|f_k\right|$. Set $C_0 = 0$ to give zero mean to the noise sequence.
  3. Set the phases of the spectral coefficients to random values, ensure symmetry if you want real valued noise.
    for $k=0..k_{\max}$
    $~~~\phi_k=\left(~\operatorname{rand\in[0,2\pi)}~\right)$,
    $~~~C_k = C_ke^{i\varphi_k}$,
    $~~~C_{-k} = C_{-k}e^{-i\varphi_k}$
  4. Take an inverse Fourier transform of the spectral coeficients to get your noise sequence, $\{y_i\} = \operatorname{ifft}(\{C_k\})$

If you want higher or lower noise power you can simply multiply the resulting sequence by a constant (or multiply all the $C_k$ values by a constant, it's the same effect).

By definition, noise generated in this way has exactly the specified spectrum since you are creating the noise by specifying its spectrum.

Improve this answer
$\endgroup$
7
  • $\begingroup$ Shouldn't $f_k=kT/(2\pi)$? $\endgroup$ – Bill Barth Feb 25 '15 at 4:17
  • $\begingroup$ Fixed, good catch @BillBarth. $\endgroup$ – Doug Lipinski Feb 25 '15 at 12:11
  • $\begingroup$ Timmer says to use gaussian complex numbers in spectral space and then scaling them, which I think is to get gaussian noise in time domain? adsabs.harvard.edu/abs/1995A%26A...300..707T $\endgroup$ – endolith Nov 11 '15 at 4:33
  • 3
    $\begingroup$ Suppose I measure a sequence of time samples from a noise process with spectrum $S(\omega)$. I then compute the discrete Fourier transform (DFT) of that time sequence and take the modulus square. The result wil not be $S(\omega)$. Due to the stochastic nature of the process, each DFT coefficient has not only a random phase but also a random amplitude. Therefore, forming frequency space points with deterministic amplitude and random phase does not faithfully reproduce the time domain statistics of a realization of the noise process. $\endgroup$ – DanielSank Jan 28 '16 at 2:22
  • $\begingroup$ @DanielSank Good point. Of course if you take enough samples the sample spectrum will converge to the true spectrum. You could add randomness to the amplitudes in some way, but it would need to depend on the number of samples and probably other things as well. I don't know the right form of noise to add to the amplitudes, but if you do then please post another answer. I think that would be very interesting. $\endgroup$ – Doug Lipinski Jan 28 '16 at 14:12
1
$\begingroup$

You can also generate noise spectrums as the solutions to stochastic differential equations. This gives an iterative method for generating the noise timestep by timestep as opposed to solving for the whole spectrum first (i.e. you just solve for one more step of the SDE, and this is your next noise term). This paper gives a good treatment of some SDEs for your noise type.

Improve this answer
$\endgroup$
1
$\begingroup$

Here's a very simple way to create pink noise in C#, which just sums lots of waves (spaced logarithmically apart) together! The functions outputs a double array with values from -1 to 1. This represents the lowest and highest points in the waveform respectively.

The 'quality' parameter represents the number of waves produced to make the sound. I find 5000 waves is just about the threshold where I can't detect any noticeable improvement with higher values, but to be on the safe side, you could (optionally) increase this to about 10,000 or higher. Also, according to Wikipedia, 20 hertz is around the lower limit of human perception in terms of what we can hear, but you can change this too if you want.

Note the sound gets quieter with a higher quality value due to technical reasons, so you may (optionally) want to adjust the volume via the volumeAdjust parameter.

public double[] createPinkNoise(double seconds, int quality=5000, double lowestFrequency=20, double highestFrequency = 44100, double volumeAdjust=1.0)
{
    long samples = (long)(44100 * seconds);
    double[] d = new double[samples];
    double[] offsets = new double[samples];
    double lowestWavelength = highestFrequency / lowestFrequency;
    Random r = new Random();
    for (int j = 0; j < quality; j++)
    {
        double wavelength = Math.Pow(lowestWavelength, (j * 1.0) / quality)  * 44100 / highestFrequency;
        double offset = r.NextDouble() * Math.PI*2;     // Important offset is needed, as otherwise all the waves will be almost in phase, and this will ruin the effect!
        for (int i = 0; i < samples; i++)
        {
            d[i] += Math.Cos(i * Math.PI * 2 / wavelength + offset) / quality * volumeAdjust;
        }
    }
    return d;
}
Improve this answer
$\endgroup$
0
$\begingroup$

I think I found a good solution:

import numpy as np
import math
import cmath
from scipy.io import wavfile


def pink_noise(f_ref, f_min, f_max, length, f_sample):
    aliasfil_len = 10000
    fil_Time = aliasfil_len * 1/f_sample
    L = length + 2 * fil_Time
    f_low = 1 / L
    f_high = f_sample
    T = f_low * 2 * np.pi
    k_max = int(f_high / f_low / 2) + 1
    print(k_max)

    # Create frequencies
    f = np.array([(k * T)/(2 * np.pi) for k in range(0, k_max)])

    # Create 1/f noise amplitude in band
    C = np.array([(1 / f[k] if (f[k] >= f_min and f[k] <= f_max) else 0)
                  for k in range(0, k_max)], dtype='complex')
    C[0] = 0
    # Create random phase in band
    Phase = np.array([(np.random.uniform(0, np.pi)
                       if (f[k] >= f_min and f[k] <= f_max)
                       else 0)
                      for k in range(0, k_max)])

    Clist_neg = list()
    Clist_pos = list()
    for k in range(-k_max + 1, -1):
        Clist_neg.append(C[-k] * cmath.exp(-1j * Phase[-k]))
    for k in range(0, k_max):
        Clist_pos.append(C[k]  * cmath.exp( 1j * Phase[k] ))

    CC = np.array(Clist_pos + Clist_neg, dtype='complex')

    # Scale to max amplitude
    maxampl = max(abs(CC))
    CC /= maxampl

    tsig = np.fft.ifft(CC)
    sig = np.real(np.sign(tsig)) * np.abs(tsig)

    # Filter aliassing
    sig = sig[aliasfil_len:-aliasfil_len]

    # clip to maximum signal and
    # correct for amplitude at reference frequency
    if f_ref > ((f_max + f_min) / 2):
        print("WARNING: f_ref ({} Hz) should be smaller or equal to the mid "
              "between {} Hz and {} Hz "
              "to prevent clipping.\n"
              "f_ref changed to {} Hz"
              .format(f_ref,
                      f_min,
                      f_max,
                      ((f_max + f_min) / 2)))
        f_ref = ((f_max + f_min) / 2)
    maxampl = max(np.abs(sig))
    sig = sig / maxampl * f_ref / ((f_max + f_min) / 2)

    halfway = int(len(sig) / 2)
    # sign invert second part for a good connection,
    # it is the mirror of the first half
    sig2nd = -1 * sig[halfway:]
    sigc = np.concatenate((sig[0:halfway], sig2nd))
    # average middle point, but second point sign inverted
    sigc[halfway] = (sig[halfway-1] - sig[halfway+1])/2

    return(sigc)


length = 10.0  # seconds
f_sample = 22050  # Hz

f_ref = 1000  # Hz, The frequency for max amplitude

f_min = 1000  # Hz
f_max = 2000  # hz

sig = pink_noise(f_ref, f_min, f_max, length, f_sample)



p = 0
for x in sig:
    p += abs(x)
rms = math.sqrt(p/len(sig))
print("rms = {:f}, {:5.1f} dB".format(rms, 20 * math.log10(rms)))

print("Length of time signal: {} samples".format(len(sig)))
print("Time signal: ", sig)
x = sig * (2**15 - 1)

wavfile.write("test.wav", f_sample, np.array(x, dtype=np.int16))
Improve this answer
$\endgroup$
1

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.