I learnt from here that the following filter can turn white noise into pink (-3dB/ octave).

  b0 = 0.99765 * b0 + white * 0.0990460;
  b1 = 0.96300 * b1 + white * 0.2965164;
  b2 = 0.57000 * b2 + white * 1.0526913;
  pink = b0 + b1 + b2 + white * 0.1848;

How would I turn white noise into blue/ azure noise (+3dB/ octave)? Is it possible to do some simple 'inversion' on this filter or can you please give me a tip as to how I can calculate the filter coefficients.


I probably should have mentioned that I need a C++ implementation but Robert's 2nd answer has me close.

Right now I do this in the header:

float *state = nullptr;

The in the implementation file:

state = new float[0.0]; in the constructor then inside the actual loop I take my white noise and do this:

float first = first_order_filter(whiteNoise, 0.99572754, 0.98443604, state);
float second = first_order_filter(first, 0.94790649, 0.83392334, state);
float third = first_order_filter(second, 0.53567505, 0.07568359, state);
out1 = third;

I am expecting pink noise as I have not yet swapped the Poles and Zeros (if I had then I would expect blue/ azure noise) but currently I get what looks and sounds like white noise.


3 Answers 3

// this processes one sample

float first_order_filter(float input, float pole, float zero, float *state)
    float new_state = input + pole*(*state);
    float output = new_state - zero*(*state);
    *state = new_state;
    return output;
  • $\begingroup$ I appreciate the help and answers but do you know what the state should be initialled with? $\endgroup$ Commented Oct 20, 2020 at 20:03
  • $\begingroup$ states should be initialized to zero. but, as long as the filters are stable ( -1 < pole < 1), then whatever initial state it is, that information will eventually die out and be forgotten by the filter. $\endgroup$ Commented Oct 20, 2020 at 23:06

You need to invert the filter, i.e. flip the poles and zeros.

The implementation that you reference is fairly awkward and will require a decent amount of math work to invert: you need to write the Z-transform for each first order section, add all the fractions into a single fraction and calculate the zeros of numerator polynomial.

An easier way would be to use this implementation https://ccrma.stanford.edu/~jos/sasp/Example_Synthesis_1_F_Noise.html and simply swap $A$ and $B$

  • $\begingroup$ That's great - I'll try that right now... $\endgroup$ Commented Oct 20, 2020 at 17:39
  • $\begingroup$ i don't think that the poles and zeros are apparent with the Kellett pink filter. but they are explicit with mine. $\endgroup$ Commented Oct 20, 2020 at 18:10
  • $\begingroup$ @robertbristow-johnson Thank you. I need a C++ implementation so I suppose I can either convert the Matlab code and swap A and B or use your poles and zeros swapped. $\endgroup$ Commented Oct 20, 2020 at 18:25
  • $\begingroup$ it's not hard to write a first-order digital filter and then cascade three of them. want me to do it? $\endgroup$ Commented Oct 20, 2020 at 18:31
  • $\begingroup$ I'll send you a Christmas card if you would! $\endgroup$ Commented Oct 20, 2020 at 18:33

here is my quarter-century old MATLAB code for a pinking filter:

% pinking filter analysis program
% from robert bristow-johnson
% please give to your local right-wing fringe group

 z = [0.98443604; 0.83392334; 0.07568359];  % the zeros   as shown in http://www.firstpr.com.au/dsp/pink-noise/
 p = [0.99572754; 0.94790649; 0.53567505];  % the poles

% z = [0.98444;    0.83392;    0.07578   ]; % the zeros   as shown in the Orfanidis book
% p = [0.99575;    0.94791;    0.53565   ]; % the poles            (these look *slightly* better)

k = 0.9297;                         % gain scaler to center error about 0 dB

[b, a] = zp2tf(z, p, k);            % convert to unfactored transfer function coefs

n = 16384;                          % number of freq points to plot

[H, w] = freqz(b, a, n);            % get frequency response

semilogx(w(2:n)/pi, 20*log10(abs(H(2:n)) + 1.0e-11));                               % plot dB gain vs. log(w)

semilogx(w(2:n)/pi, 20*log10(abs(H(2:n)) + 1.0e-11) + 20*log10(sqrt(w(2:n)) + 1.0e-11));        % plot dB error vs. log(w)

axis([0.0005 1 -0.35 0.35]);            % close up view of the error

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.