Making grey audio noise using matlab / octave

Question

I can create pink,brown,blue audio noise with the following snippet code by selecting a different invfnorm variable to use below, but how can I create grey noise?

%https://en.wikipedia.org/wiki/Colors_of_noise

mean_amp=mean(yamp_orig.^2); %get mean of all freq amplitudes
amt_of_freq=size(xfreq_orig,1); %number of freq

%invfnorm=1./[1:amt_of_freq]; % 1/f creates pink noise
%invfnorm=[1:amt_of_freq]; % f creates blue noise
invfnorm=1./[1:amt_of_freq].^2; % 1/f^2 creates brown noise

amp_1f_new=sqrt(mean_amp*invfnorm/sum(invfnorm))(:); %new noise amplitudes to use

In the link https://en.wikipedia.org/wiki/Colors_of_noise they give the formulas for pink,brown,blue audio noise but for the grey noise they just say it's "an inverted A-weighting curve" not showing the formula, I just need the formula. see spectrum below.

The website I got the information for these are located here grey noise

Ps: I'm using Octave 4.2.2 which is similar to matlab

Answer 1

As you mentioned, grey noise is created by applying an inverse a-weighting curve.

The following snippet is a Matlab example (thanks to W. Owen Brimijoin) for generating grey noise:

%values from the ISO 66-phon Equal-loudness contour (adjusted for
        %optimal spline interpolation):
        freqs = [1 5 15 36 75 138 235 376 572 835 1181 1500 2183 2874 3718 ...
            4800 5946 7377 9051 10996 13239 15808 18735 22050].*sample_rate/44100;
        dB_vals = [61 61 56 40 25 17 11 7 5 4 6 8 3 1 1 4 9 14 17 16 10 5 2 1];

        %create level vector for use in inverse Fourier transform:
        freq = linspace(0,sample_rate/2,floor(num_samples/2));
        spl = spline(freqs,dB_vals,freq); %upsample 
        levels = [spl,fliplr(spl)]; %reflect vector
        levels = 10.^(levels'./10); %change to power
        levels(levels==inf) = 0; %remove infinite values
        phase_vals = rand(length(levels),1); %generate random phase vals
        %now apply an inverse fft:
        wave = real(ifft(sqrt(levels).*(cos(2*pi*phase_vals)+1i*sin(2*pi*phase_vals))));
        wave = wave./max(abs(wave));

Where levels = the inversed a-weighting array.

This example creates a grey-noise oscillator and applies the filter in the frequency domain to the signal, before turning the signal back to the time domain.

As you mentioned you said you just needed the formula so maybe this specific line helps the most (once you have already calculated the inversed a-weighting curve) :

wave = real(ifft(sqrt(levels).*(cos(2*pi*phase_vals)+1i*sin(2*pi*phase_vals))));

As stated at the start of the code, this is using the 66 phon curve so you may need to adjust your dB_vals array if you wanted to use a different level phon.

I found this function quite useful for calculating various phon curves for a-weighting, if necessary.