A Study in Pink

Programming and Development Nyquist
1

Steve asked me lately if I wouldn’t know an easier method to produce pink noise (See his method in the “Add-noise”-plug-in).
Why do we want a method for this kind of noise, isn’t it not already available in the generate menu?
Yes, it is, but Nyquist has no access to this C++ code.
The only thing that is available is white noise.
Coloured noise is white noise of which the amplitude decreases the more, the higher the frequency gets. Thus filtering with a lowpass is the logical choice.
Something like

(LP (NOISE) 20)

Ok, we have a sound with attenuated higher frequencies. Unfortunately, the low pass with the lowest order has a roll off of 6 dB per octave, which produces 1/f^2 noise, so called “brown” noise.
Pink noise is “1/f”, this means that the roll off or transition band should be damped by 3 dB per octave.
Although no filter with that roll off is available, here’s a bi-quad implementation that comes fairly close:

;; original first 3 coeffs b0-2 were 5 times smaller
;; Fine for measurement with sine, 
;; but 14 dB too small resulting noise
(defun pink-lp (snd)
  (biquad-m snd
    0.247876311 -0.315279067 0.074161016
    1.0 -1.80116083982126 0.80257737639225))

(pink-lp (noise) )

(I’ve found the original code in one of Edgar’s plug-ins)
The measurement with sine waves returns fairly good damping values of 3 dB as an average.
But if the LP is applied as such on white noise, the output is about 14 dB too low, don’t ask me why. I therefore multiplied the first 3 coeffs by 5. Let me know if the output is too high on your system (all is possible with 3 different OS’ and various architectures).

2

That’s pretty close.
The coefficients used by Edgar appear to have originally come from here: http://www.musicdsp.org/

3

Yes, that’s the place.
There’s also a first order filter, but it isn’t accurate enough, thus I haven’t listed it.
a cool thing is also the highpass-version of the above mentioned filter. Thus you can create lowpass filters with 9/15/21 etc. dB roll-off. That’s especially useful for control signals of any kind.

4

Here’s a function that generates pink noise without filtering. The algorithm is the same as is used in the Audacity pink noise generator.
Unfortunately it is rather slow, but the quality is good.

;;; Generate pink noise
(defun pink (&optional (dur 1.0)(amp 1.0))
  (setf pink-class (send class :new '(dur)))
  (send pink-class :answer :next '()
  '((setq white (1- (* 2 (rrandom))))
    (psetq buf0 (+ (* 0.997 buf0)(* 0.029591 white))
           buf1 (+ (* 0.985 buf1)(* 0.032534 white))
           buf2 (+ (* 0.950 buf2)(* 0.048056 white))
           buf3 (+ (* 0.850 buf3)(* 0.090579 white))
           buf4 (+ (* 0.620 buf4)(* 0.108990 white))
           buf5 (+ (* 0.250 buf5)(* 0.255784 white))
           dur (1- dur))
      (if (>= dur 0)
          (* 0.55 (+ buf0 buf1 buf2 buf3 buf4 buf5))
          nil)))
  (send pink-class :answer :isnew '(p1)
    '((setf dur p1)))
  (defun pink-mono (dur amp)
    (psetq buf0 0 buf1 0 buf2 0
           buf3 0 buf4 0 buf5 0)
    (let (obj)
      (setf obj (send pink-class :new dur))
      (mult amp (snd-fromobject 0 *sound-srate* obj))))
  (cond
    ((or (soundp s)(arrayp s))
       (setq dur (* dur len)))
    (t (setq dur (* dur *sound-srate*))
       (setq len dur)))
  (multichan-expand #'pink-mono dur amp))


(setq duration 1)
(setq amplitude 0.5)
(pink duration amplitude)
5

It is actually nothing but white noise that is filtered too (see the variable “white”).
But instead of one biquad of II.-order, there are 6 parallel I.-order Biquads used.

The fact that we can use the biquad functions in this case too, allows us to increase the speed by maybe 50 times. However, the beautifully coded object-oriented snippet from above is replaced by a very profane one. The only thing that I’ve changed in comparison to the original C++ code is the normalisation value (0.55), which is therin described as “experimental”.
A lot of tests (about 100000 for each duration value of 0.0001, 0.001 … 10000.0 seconds ) have shown that the peak depends on the duration length. I’ve included a routine that interpolates the different scale-factors. Although the normalisation should be better than the noise in the generate menu, it is still possible that the amplitude exceeds 0 dB, we’re dealing with random noise after all.
Furthermore, I can’t say how much the round-off errors influence the results.

(set-sound-srate (setf *sr* *sound-srate*))
(defun snd-pink (dur &optional (srate *sound-srate*))
  (let*
        ((white (snd-white 0 srate dur))  
         (norm '((10000 0.335) (1000 0.35) (100 0.36)
                 (10 0.39) (1 0.457) (0.1 0.532)
                 (0.01 0.680) (0.001 1.008) (0.0001 1.906)))
         (pair1 (assoc dur  norm :test '>=))
         (pair2 (assoc dur (reverse norm) :test '<=))
         (scale-factor (if (equal pair1 pair2) (second pair1) 
            (+ (second pair1) (* (- dur  (first pair1))
                       (/ (- (second pair2)  (second pair1)) 
                       (- (first pair2)(first pair1)))))))
         (coeffs '((0.997 0.029591) (0.985 0.032534)
                   (0.950 0.048056) (0.850 0.090579)
                   (0.620 0.108990) (0.250 0.255784))))
        (mult scale-factor
              (simrep (i 6)
                 (biquad white (cadr (nth i coeffs)) 0 0 
                               1 (car (nth i coeffs)) 0)))))
;; test noise
(snd-pink 20)
6

What I meant was that it does not use any of the Nyquist filter functions.
I thought that it should be possible to rewrite in terms of biquad filters, but I was unsure how to do that. Your code makes that transposition nicely, and yes it’s a lot faster.

So long as very long tracks are not required, the output could simply be normalized:

(defun snd-pink (dur amp)
  (setf s
    (let ((white (abs-env (noise dur))))
      (sim
        (biquad white 0.029591 0 0 1 0.997 0)
        (biquad white 0.032534 0 0 1 0.985 0)
        (biquad white 0.048056 0 0 1 0.950 0)
        (biquad white 0.090579 0 0 1 0.850 0)
        (biquad white 0.108990 0 0 1 0.620 0)
        (biquad white 0.255784 0 0 1 0.250 0))))
  (let ((norm (/ (peak s ny:all))))
    (mult amp norm s)))

(defun pink (&optional (dur 1) (amp 1))
  (multichan-expand #'snd-pink dur amp))

(pink 30)
7

Filter design is always tricky. Even if one knows the coefficients, you still do not know after which convention they have to be passed to biquad or biquad-m.
It took me a while to figure it out although it now appears so simple.
How is the quality? Do you have an idea how to test the result? There are many algorithms that use the same approach with slightly different numbers and less or more filters. But I have no clue how to compare them.

Ideally, the above algorithm should be able to produce all coloured noise (1/F^alpha) but I guess it is much too hard to re-calculate all coefficients.
Can you think of other noise types that we could implement?
I am especially interested in chaotic noise (such as the lorenz system). But for those, the sample by sample method may be more appropriate.

8

Seems pretty good.
One measure is the spectrum, which can be measured with “Plot Spectrum”. It should also sound smooth and even. It’s quite interesting that if you clip pink noise, even quite severely, there is little change to the spectrum, but the sound changes and becomes kind of “lumpy”.

Applying a dB per decade FFT filter to white noise should give pretty good pink noise, and the filter can be created fairly easily using the Equalization effect - the sound of this can then be compare with our generated pink noise.

Pink noise is probably the most useful colour.

Brown noise might be useful. The Audacity code for brown noise is pretty simple and should be easy to implement in Nyquist:

   case 2: // brown
       // fc=100 Hz,
       // y[n]=a0*x[n] + b1*y[n-1];
       white=a0=b1=fc=y=0;
       fc=100; //fs=44100;
       b1=exp(-2*M_PI*fc/fs);
       a0=1.0f-b1;

       // 6.2f is an experimental normalization factor: thanks to Martyn
       amplitude *= 6.2f;
       for(i=0; i<len; i++){
         white=(rand() / div) - 1.0f;
         y = (a0 * white + b1 * y);
         buffer[i] = amplitude * y;
       };
9

It seems that brown noise is implemented as one I.-order biquad filter. The cut-off frequency is fairly high though.
In general, I would simply use the normal LP function with a cut-off at about 20 Hz.
The great question is which Lowpass (with 6 dB per octave) is suited the best.
Normally, it is most desireable to have a flat freq response in the pass band (=Butterworth), but in the case of brown noise, it is more important that the transition band is stable (i.e. has no ripples). That’s why I search for a method to test the filter.
Lowpasses tend to be very instable when the cut-off frequency is near the lower limit. Furthermore, the filters in Nyquist seem all to be designed for sampling rates of 44100 Hz, thus an adaption of the above source code would be logical.
In principle, it would be very easy to get a desired roll-off damping by simply mixing unprocessed noise with the same noise that is filtered (e.g. 1/2 times 0 dB + -6 dB = -3 dB). However, that’s only possible if the phase response is ideal.
After all, we do not have to be too much concerned about perfect results since critical mathematical algorithms can’t be achieved within Nyquist. The floating point precision is simply too bad and messy (sometimes it is float sometimes double).
For the purpose of sound synthesis, our achieved results seem to be fairly good.
I wonder how “white” our white noise itself is, because all filtered noise depends on the source material.

10

Apropos white noise. Here are two procedures that generate Gauss-distributed noise, i.e. most values appear around 0 (top of the bell curve).
The noise is less hissing than white noise. I’ve come up with two different approaches but the results seem to be essentially the same.

(defun snd-gauss2 (dur)
  (let* 
        ((gauss (simrep (i 12) (snd-white 0 *sound-srate* dur)))
         (norm-scale (peak gauss ny:all)))
        (scale (/ norm-scale)  gauss)))

(defun snd-gauss (dur)
  (let* 
        ((gauss (snd-avg
                         (mult 1.414 
                               (snd-white 0 (* 12 *sound-srate*) dur))
                         12 12 op-average))
         (norm-scale (peak gauss ny:all)))
        (scale (/ norm-scale)  gauss)))
(setf white (scale 0.1 (abs-env (noise 2))))
(seq white 
     (mult 0.1 (snd-gauss 2))
     (mult 0.1 (snd-gauss2 2)))

It is said that random variables are gaussian distributed when they are at least 12 times added up, that’s the principle behind these functions. I can’t say which one is better or faster, both took about the same amount of time for 1000 seconds (~14 s).

11

An alternative approach that gives rather nice brown noise with 6dB/octave all the way down to 20 Hz.

(defun brown (&optional (dur 1)(amp 1))
  (multichan-expand #'snd-brown dur amp))

(defun snd-brown (dur amp)
  (let* ((white (abs-env (noise dur)))
         (brown (highpass8 (integrate white) 15))
         (gain (/ (peak brown ny:all))))
    (mult amp gain brown)))

(brown 20 0.5)
12

There does not seem to be much difference in the sound if the rms level is the same.
This audio sample has been created in three parts, your two Gaussian code snippets and the Audacity white noise generator,
The peak amplitude of the Gaussian noise is about 7.7 dB higher than the Audacity white noise.

13

There shouldn’t be much difference, it is still white noise after all and my bad speakers may exhibit a total different spectrum from yours.
It gets more prominent after filtering with a resonance filter or after a hard clipping.
In principle, the random values should be analysed statistically to see how they are distributed before/after (although the RMS may be a good indicater as it is).
Besides, interesting brown filter snippet.

14

I think I’m getting the hang of this biquad stuff now :smiley:

So here is a really inefficient (slow) way to generate bandwidth limited brown noise:

(setf sarray (make-array 44100))
(setq hz 20 q 0.707)
(setq w (* 2.0 Pi (/ hz *sound-srate*))
      cw (cos w)
      sw (sin w)
      alpha (* sw (sinh (/ 0.5 q)))
      a0 (+ 1.0 alpha)
      a1 (* -2.0 cw)
      a2 (- 1.0 alpha)
      b1 (- -1.0 cw)
      b0 (* -0.5 b1)
      b2 b0)

(defun zrand ()
  (* 2 (- (rrandom) 0.5)))

(defun brown (amp)
  (setq maxval 0)
  (do ((i 0 (1+ i)) (in (zrand)(+ in (zrand)))
       (x1 0) (x2 0) (y1 0) (y2 0) (in 0) (out 0))
      ((= i (length sarray)))
    (setq out (/ (- (+ (* b0 in) (* b1 x1) (* b2 x2))(* a1 y1) (* a2 y2)) a0))
    (setq x2 x1)
    (setq x1 in)
    (setq y2 y1)
    (setq y1 out)
    (when (> (abs out) maxval)(setq maxval (abs out)))
    (setf (aref sarray i) out))
  (mult (/ maxval)(snd-from-array 0 *sound-srate* sarray)))

(brown 0.8)

and the point of writing such inefficient code: it is only slow because it is looping in LISP. If written in C++ it would be very much faster, so it could be used to improve the Audacity brown noise generator.

15

Can’t you use the normal Biquad function with the calculated values?
Or are there some coeffs updated during execution? Because, that’s the only negative thing about Biquad, it can only work with constants (in contrast to the built-in LP and HP functions, where at least the cut-off can be changed).
I feel a slight headache coming up when I am looking at your code… :smiley:

16

The coefficients change depending on the sample rate. By calculating them at run time the filter frequency can be kept the same.

At least I didn’t write it object-oriented :smiley:

17

You’ve forgotten to amplify the brown noise by the “amp” factor.
I guess the last line of the function should be

  (mult (/ amp maxval)(snd-from-array 0 *sound-srate* sarray)))

I still think that the whole thing could be written in terms of a biquad function.
The input sample is obviously continously summed with a random number. (besides, “in” is two times initialized in the do loop).
Therefore, the input sound for “Biquad” should be integrated white noise (multiplied by
sound-srate). Maxval would be the peak of this sound. Do you agree so far?

18

Not “forgotten”, I “didn’t bother”. It’s not stereo either, but it works well enough as a proof of concept.

Anyway, I think you’ll like this. It’s based on Paul Kellet’s “instrumentation grade” pink noise.

(defun pink (dur amp)
  (mult amp (db-to-linear (- (+ 15 (power dur 0.1))))
    (let ((white (abs-env (noise dur))))
      (sim
        (biquad white 0.0555179 0 0 1 0.99886 0)
        (biquad white 0.0750759 0 0 1 0.99332 0)
        (biquad white 0.1538520 0 0 1 0.96900 0)
        (biquad white 0.3104856 0 0 1 0.86650 0)
        (biquad white 0.5329522 0 0 1 0.55000 0)
        (mult -1 (biquad white 0.0168980 0 0 1 -0.7616 0))))))

(pink 10 1)
19

Ok, I’ll take a look.

Here’s your brown filter in comparison to the one with biquad. I’ve limited the sound to 10 places in order to facilitate the comparing of the two results. As I mentioned before, the in-sound could be integrated white noise. (full functional code at the end of the post)

(setf sarray (make-array 10))
(setf zrand (vector 0.5 0.2 0.3 -0.1 -0.7 -1.2 -0.6 0.1 0.8 1.3))
(setq hz 20 q 0.707)
(setq w (* 2.0 Pi (/ hz *sound-srate*))
      cw (cos w)
      sw (sin w)
      alpha (* sw (sinh (/ 0.5 q)))
      a0 (+ 1.0 alpha)
      a1 (* -2.0 cw)
      a2 (- 1.0 alpha)
      b1 (- -1.0 cw)
      b0 (* -0.5 b1)
      b2 b0)

(defun brown (amp)
  (setq maxval 0)
  (do ((i 0 (1+ i)) 
       (x1 0) (x2 0) (y1 0) (y2 0) (out 0))
      ((= i (length sarray)))
    (setf in (aref zrand i))
    (setq out (/ (- (+ (* b0 in) (* b1 x1) (* b2 x2))(* a1 y1) (* a2 y2)) a0))
    (setq x2 x1)
    (setq x1 in)
    (setq y2 y1)
    (setq y1 out)
    (when (> (abs out) maxval)(setq maxval (abs out)))
    (setf (aref sarray i) out))
  (mult (/ amp maxval)(snd-from-array 0 *sound-srate* sarray)))

(setf amp 0.8)
(setf snd-in (snd-from-array 0 44100 zrand) )
(setf out (biquad-m  snd-in b0 b1 b2 a0 a1 a2))
(setf maxval (peak out ny:all))
(snd-display (mult (/ amp maxval) out))
(snd-display (brown amp))

I wonder if it is enough to normalise the result by only one positive or negative value, or if it is necessary to remove the DC?

Ps. Here’s my version:

(defun snd-brown (&optional (amp 1.0) (dur 1.0))
  (let*
        ((hz 20)
         (q (/ (sqrt 2.0) 2.0))
         (w (* 2.0 Pi (/ hz *sound-srate*)))
         (cw (cos w))
         (sw (sin w))
         (alpha (* sw (sinh (/ 0.5 q))))
         (a0 (+ 1.0 alpha)) (a1 (* -2.0 cw)) (a2 (- 1.0 alpha))
         (b1 (- -1.0 cw)) (b0 (* -0.5 b1)) (b2 b0)
         (white (mult *sound-srate* 
                      (integrate (abs-env (noise dur)))))
         (out (biquad-m  white b0 b1 b2 a0 a1 a2))
         (maxval (peak out ny:all)))
        (mult (/ amp maxval) out)))
(defun s-brown (&optional (amp 1.0) (dur 1.0)) 
  (multichan-expand 'snd-brown amp dur))
(s-brown 0.8 2.0)
20

Thinking about brown noise and the impossibility to predict the peak output level, I was wondering how much damage would be done by constraining the sample values.
Clipping brown noise to a specified range has quite a bad effect on the sound, but what about if peaks that drift beyond the target level are simply “folded over” back into range?

Brown noise can be defined in terms of adding a random number to the previous sample value. On average the positive and negative random offsets are likely to balance out reasonably well and the net offset would be expected to be fairly small, but if the offsets are truly random, then there is a chance that the net offset could be huge. This behaviour can be ameliorated by removing DC offset, but there is still a chance that there could be a huge positive going drift, followed by an equally huge negative drift.

Another approach to keeping brown noise within a reasonable range is to filter out low frequencies. For audio use it would seem reasonable to roll-off frequencies below 20 Hz (as the bottom end of the audio spectrum). This works pretty well, but it assumes that the user does not want low frequency brown noise, but a user may want low frequencies, for example to use as a control signal or may not be dealing with audio but some other type of signal.

The solution that I’m considering is to “folder over” any peaks that are beyond 0 dB, so for example if we have sample values:
0.9 1.0 1.1 1.2 1.0 0.9
then they become:
0.9 1.0 0.9 0.8 1.0 0.9

Inverted brown noise is still brown noise, so as long as we don’t mess with the randomness too often it should make little difference to the quality of the noise.
Of course there is still a chance that the amplitude could be lower than expected, but when playing with random events, that’s the chance you take.

A short snippet to illustrate “folding over”:

(let ((top (sum 2 (mult -1 (s-max s 1))))
      (bottom (sum -2 (mult -1 (s-min s -1)))))
  (sim (clip s 1) top bottom))

This is not ideal as it messes with the randomness a lot more than necessary, but in practice it seems to work reasonably well, and certainly sound better than clipping.

One last point: If brown noise decreases in power by 6 dB per octave, then for ideal brown noise, what is the DC offset?