Pink noise is very useful in computer music and audio engineering, and I have always wanted an accurate pink noise object in the Pure Data environment. Generating accurate pink noise is not as simple as it might seem, so I spent some time designing a so-called pinking filter, and generalized it for any audio sample rate. This page describes the filter I ended up with, shows how its parameters are calculated, and displays the Pd object source code at the bottom of the page.

Download
pink~.zip — Pd object (os x), help patch, source code, and matlab parameter calculator

Algorithm
The pinking filter creates pink noise by filtering white noise to create the desired pink energy spectrum. The filter's frequency response is -3dB per octave starting at approximately four Hertz and continuing to just below the Nyquist frequency. The pinking filter consists of several one-pole low-pass filters, where each low-pass is spaced two octaves from its neighbors and filter gains are attenuated in 6dB steps. This configuration results in a nearly linear -3dB per octave overall frequency response when the low-pass filters are summed.

This algorithm is designed to work with any sample rate up to 192kHz. Filter coefficients are calculated by first determining the maximum cutoff frequency resulting in a filter coefficient less than one, described by the following formula where f is frequency and sr is sample rate. The - 1 ensures that the coefficient resulting from f is just below 1 to prevent the first filter from exploding.
f = sr / 2π - 1
Following this, a series of low-pass coefficents are calculated by successively dividing frequency by four while it remains larger than one. This operation is expressed by the following while loop, where a is an array of coefficients and i is the array index starting at zero.
while( f > 1 )
{
    a[i] = 2π * f / sr
    f    = f / 4
    i++
}
The low-pass filter gains are calculated starting at 0dB for the lowest frequency filter, and continue downward in 6dB steps for each successive filter. An additional output gain factor is calculated so that the summed low-pass filters result in an overall frequency response that starts at 0dB-FS before proceeding downward. This process is described by the following while loop where g is an array of gain values, dB is the decibel gain value starting at zero, s is the sum of gain values starting at zero, og is the output gain factor, and i is the array index starting at the number of filters.
while( i-- > 0 )
{
    g[i]  = pow( 10, dB / 20 )
    s    += g[i];
    dB   -= 6.0;
}

og = pow( 10, ( ( -20 * log10( s ) ) / 20 )
Finally, pink noise is created by filtering white noise through the low-pass filters. The filters are implemented within a loop so that a variable number of filters may be used depending on sample rate, and this loop runs once for each sample of output. This implementation is expressed by the following statements where w is the white noise sample between -1 and 1, p is the pink noise sample initialized to zero, i is the array index starting at -1, f is the number of filters, y is the array of filtered samples, a is the array of filter coefficients, g is the array of filter gains, o is the output sample, and og is the output gain.
while( ++i < f )
{
    y[i] = a[i] * w + ( 1 - a[i] ) * y[i];
    p   += y[i] * g[i];
}

o = og * p;
The graph below shows the pinking filter's frequency response starting at 0dB, as well as the individual low-pass filters' frequency responses below that. Sample rate is 44.1kHz.


Cutoff Frequencies, Coefficients, and Gains
Below are lists of filter pararameters for common audio sample rates. Parameters for other sample rates may be calculated with the pink.m matlab script in pink~.zip.

Sample Rate : 8000
Output Gain : 0.506845730

Freq            Coeff           Gain
-------------------------------------------
1.242421430     0.000975796     1.000000000
4.969685722     0.003903182     0.501187234
19.878742886    0.015612728     0.251188643
79.514971546    0.062450913     0.125892541
318.059886184   0.249803650     0.063095734
1272.239544735  0.999214602     0.031622777
    
Sample Rate : 11025
Output Gain : 0.506845730

Freq            Coeff           Gain
-------------------------------------------
1.712581296     0.000976006     1.000000000
6.850325186     0.003904024     0.501187234
27.401300744    0.015616095     0.251188643
109.605202974   0.062464381     0.125892541
438.420811897   0.249857524     0.063095734
1753.683247588  0.999430097     0.031622777
Sample Rate : 16000
Output Gain : 0.506845730

Freq            Coeff           Gain
-------------------------------------------
2.485819423     0.000976179     1.000000000
9.943277693     0.003904716     0.501187234
39.773110773    0.015618864     0.251188643
159.092443092   0.062475456     0.125892541
636.369772368   0.249901825     0.063095734
2545.479089470  0.999607301     0.031622777
    
Sample Rate : 22050
Output Gain : 0.506845730

Freq            Coeff           Gain
-------------------------------------------
3.426139155     0.000976284     1.000000000
13.704556622    0.003905137     0.501187234
54.818226487    0.015620548     0.251188643
219.272905949   0.062482191     0.125892541
877.091623794   0.249928762     0.063095734
3508.366495176  0.999715048     0.031622777
Sample Rate : 32000
Output Gain : 0.502806702

Freq            Coeff           Gain
-------------------------------------------
1.243153852     0.000244093     1.000000000
4.972615409     0.000976371     0.501187234
19.890461636    0.003905483     0.251188643
79.561846546    0.015621932     0.125892541
318.247386184   0.062487728     0.063095734
1272.989544735  0.249950913     0.031622777
5091.958178941  0.999803650     0.015848932
    
Sample Rate : 44100
Output Gain : 0.502806702

Freq            Coeff           Gain
-------------------------------------------
1.713313718     0.000244106     1.000000000
6.853254873     0.000976423     0.501187234
27.413019494    0.003905693     0.251188643
109.652077974   0.015622774     0.125892541
438.608311897   0.062491095     0.063095734
1754.433247588  0.249964381     0.031622777
7017.732990353  0.999857524     0.015848932
Sample Rate : 48000
Output Gain : 0.502806702

Freq            Coeff           Gain
-------------------------------------------
1.864852849     0.000244109     1.000000000
7.459411395     0.000976435     0.501187234
29.837645580    0.003905739     0.251188643
119.350582319   0.015622955     0.125892541
477.402329276   0.062491819     0.063095734
1909.609317103  0.249967275     0.031622777
7638.437268411  0.999869100     0.015848932
    
Sample Rate : 64000
Output Gain : 0.502806702

Freq            Coeff           Gain
-------------------------------------------
2.486551845     0.000244117     1.000000000
9.946207381     0.000976467     0.501187234
39.784829523    0.003905867     0.251188643
159.139318092   0.015623466     0.125892541
636.557272368   0.062493864     0.063095734
2546.229089470  0.249975456     0.031622777
10184.916357881 0.999901825     0.015848932
Sample Rate : 88200
Output Gain : 0.502806702

Freq            Coeff           Gain
-------------------------------------------
3.426871577     0.000244123     1.000000000
13.707486309    0.000976493     0.501187234
54.829945237    0.003905972     0.251188643
219.319780949   0.015623887     0.125892541
877.279123794   0.062495548     0.063095734
3509.116495176  0.249982191     0.031622777
14036.465980705 0.999928762     0.015848932
    
Sample Rate : 96000
Output Gain : 0.502806702

Freq            Coeff           Gain
-------------------------------------------
3.729949838     0.000244125     1.000000000
14.919799352    0.000976499     0.501187234
59.679197409    0.003905994     0.251188643
238.716789638   0.015623977     0.125892541
954.867158551   0.062495909     0.063095734
3819.468634205  0.249983638     0.031622777
15277.874536822 0.999934550     0.015848932
Sample Rate : 176400
Output Gain : 0.500806513

Freq            Coeff           Gain
-------------------------------------------
1.713496824     0.000061033     1.000000000
6.853987295     0.000244132     0.501187234
27.415949181    0.000976528     0.251188643
109.663796724   0.003906111     0.125892541
438.655186897   0.015624443     0.063095734
1754.620747588  0.062497774     0.031622777
7018.482990353  0.249991095     0.015848932
28073.931961410 0.999964381     0.007943282
    
Sample Rate : 192000
Output Gain : 0.500806513

Freq            Coeff           Gain
-------------------------------------------
1.865035954     0.000061033     1.000000000
7.460143817     0.000244133     0.501187234
29.840575267    0.000976531     0.251188643
119.362301069   0.003906122     0.125892541
477.449204276   0.015624489     0.063095734
1909.796817103  0.062497955     0.031622777
7639.187268411  0.249991819     0.015848932
30556.749073644 0.999967275     0.007943282

References
The technique described above was inspired by Paul Kellet's filter coefficients, located here:
· Filter to make pink noise from white

Larry Trammell's method is described here:
· Pink Noise Generator
· Trammell Pink Noise (C++ class)

An excellent reference page about pink noise techniques is located here:
· DSP generation of Pink (1/f) Noise

My exploratory Pd object source code for the Kellet and Trammel techniques is available here:
· pink~.ref.zip

pink~.c

Deprecated: Function create_function() is deprecated in /home/c00p/www/lib/geshi/geshi.php on line 4698
//------------------------------------------------------------------------------
//  Pink Noise Object for Pd
//
//  pink~.c
//
//  Generates pink noise using a series of summed and attenuated lowpass filters
//
//  Created by Cooper Baker on 2/4/15.
//  Copyright (c) 2015 Cooper Baker. All rights reserved.
//------------------------------------------------------------------------------


//------------------------------------------------------------------------------
// headers
//------------------------------------------------------------------------------

// main header for pd
#include "m_pd.h"

// standard library for random number generator
#include <stdlib.h>

// time header to seed random numbers
#include <time.h>

// utility header for Pd Spectral Toolkit project
#include "utility.h"

// disable compiler warnings on windows
#ifdef NT
#pragma warning( disable : 4244 )
#pragma warning( disable : 4305 )
#endif


//------------------------------------------------------------------------------
// pink_class - pointer to this object's definition
//------------------------------------------------------------------------------
static t_class* pink_class;


//------------------------------------------------------------------------------
// pink - data structure holding this object's data
//------------------------------------------------------------------------------
typedef struct pink
{
    // this object - must always be first variable in struct
    t_object object;
   
    // filter coefficients
    double a[ 8 ];

    // filter gains
    double g[ 8 ];
   
    // filter states
    double y[ 8 ];

    // filter gain
    double gain;
   
    // number of filters
    t_int filters;
   
    // random normalization coefficient
    double rand_norm;
   
} t_pink;


//------------------------------------------------------------------------------
// function prototypes
//------------------------------------------------------------------------------
static t_int* pink_perform     ( t_int* io );
static void   pink_dsp         ( t_pink* object, t_signal **sig );
static void*  pink_new         ( void );
void          pink_tilde_setup ( void );


//------------------------------------------------------------------------------
// pink_perform - the signal processing function of this object
//------------------------------------------------------------------------------
static t_int* pink_perform( t_int* io )
{
    // store variables from dsp input/output array
    t_float* out    = ( t_float* )( io[ 1 ] );
    t_int    frames = ( t_int    )( io[ 2 ] );
    t_pink*  object = ( t_pink*  )( io[ 3 ] );

    // copy pointers to filter value arrays
    double* a = object->a;
    double* g = object->g;
    double* y = object->y;
   
    // copy filter parameters
    double gain    = object->gain;
    t_int  filters = object->filters;
   
    // random number variable
    double rand_norm = object->rand_norm;

    // signal vector iterator variable
    t_int n = -1;
   
    // filter iterator variable
    t_int i;
   
    // noise variables
    double white;
    double pink;
   
    // the dsp loop
    while( ++n < frames )
    {
        // make white noise normalized between -1 and 1
        white = rand() * rand_norm - 1.0;

        // initialize
        pink  = 0.0;
        i     = -1;

        // the filter loop
        while( ++i < filters )
        {
            // filter the white noise
            y[ i ] = a[ i ] * white + ( 1.0 - a[ i ] ) * y[ i ];
           
            // apply gain and accumulate filtered noise
            pink += y[ i ] * g[ i ];
        }
       
        // apply overall gain and copy to output
        out[ n ] = pink * gain;
    }

    // return the dsp input/output array address plus one more than its size
    // to provide a pointer to the next perform function in pd's call list
    return &( io[ 4 ] );
}


//------------------------------------------------------------------------------
// pink_dsp - installs this object's dsp function in pd's callback list
//------------------------------------------------------------------------------
static void pink_dsp( t_pink* object, t_signal **sig )
{
    // get sample rate
    t_float sr = sig[ 0 ]->s_sr;

    // copy pointers to filter value arrays
    double* a = object->a;
    double* g = object->g;

    // allocate temporary variables
    t_int  i       = 0;
    t_int  db      = 0;
    double amp_sum = 0;
   
    // calculate maximum cutoff frequency so that filter coefficient < 1.0;
    double freq = ( sr / C_2_PI ) - 1.0;
   
    // calculate coefficients
    while( freq > 1 )
    {
        a[ i ] = C_2_PI * freq / sr;
        freq   = freq / 4.0;
       
        i++;
    }

    // store number of filters
    object->filters = i;
   
    // calculate gains
    while( i-- > 0 )
    {
        g[ i ]   = DbToA( db );
        amp_sum += DbToA( db );
        db      -= 6.0;
    }
   
    // calculate overall gain
    object->gain = DbToA( -AToDb( amp_sum ) );
   
    // dsp_add arguments
    //--------------------------------------------------------------------------
    // perform routine
    // number of passed parameters
    // outlet sample vector
    // sample frames to process (vector size)
    // object pointer
    dsp_add( pink_perform, 3, sig[ 0 ]->s_vec, sig[ 0 ]->s_n, object );
}


//------------------------------------------------------------------------------
// pink_new - instantiates a copy of this object in pd
//------------------------------------------------------------------------------
static void* pink_new( void )
{
    // create a pointer to this object
    t_pink* object = ( t_pink* )pd_new( pink_class );
   
    // create a new signal outlet for this object
    outlet_new( &object->object, gensym( "signal" ) );
   
    // calculate 2 over rand max
    object->rand_norm = ( 2.0 / RAND_MAX );

    // initialize filter values
    memset( object->a, 0, 8 * sizeof( double ) );
    memset( object->g, 0, 8 * sizeof( double ) );
    memset( object->y, 0, 8 * sizeof( double ) );
   
    // seed the random number generator
    time_t t;
    srand( ( unsigned )time( &t ) );
   
    return object;
}


//------------------------------------------------------------------------------
// pink_tilde_setup - describes the attributes of this object to pd so it may be properly instantiated
// (must always be named with _tilde replacing ~ in the object name)
//------------------------------------------------------------------------------
void pink_tilde_setup( void )
{
    // creates an instance of this object and describes it to pd
    pink_class = class_new( gensym( "pink~" ), ( t_newmethod )pink_new, 0, sizeof( t_pink ), CLASS_NOINLET, 0 );
   
    // installs pink_dsp so that it will be called when dsp is turned on
    class_addmethod( pink_class, ( t_method )pink_dsp, gensym( "dsp" ), 0 );
   
    // announce this object in the pd console
    Announce( "pink~: pink noise generator - v1.0" );
}


//------------------------------------------------------------------------------
// EOF
//------------------------------------------------------------------------------