# Given a -3dB/ octave filter that makes Pink Noise, how can I make a +3dB/ octave filter that makes Blue/ Azure Noise?

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.

EDIT

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.

// 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;
}

• I appreciate the help and answers but do you know what the state should be initialled with? – Dave Chambers Oct 20 '20 at 20:03
• 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. – robert bristow-johnson Oct 20 '20 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$$

• That's great - I'll try that right now... – Dave Chambers Oct 20 '20 at 17:39
• i don't think that the poles and zeros are apparent with the Kellett pink filter. but they are explicit with mine. – robert bristow-johnson Oct 20 '20 at 18:10
• @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. – Dave Chambers Oct 20 '20 at 18:25
• it's not hard to write a first-order digital filter and then cascade three of them. want me to do it? – robert bristow-johnson Oct 20 '20 at 18:31
• I'll send you a Christmas card if you would! – Dave Chambers Oct 20 '20 at 18:33

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

%
% pinking filter analysis program
%
% from robert bristow-johnson
%
%

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)
pause;

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)
pause;

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