# (Dis-)connecting a stochastic signals spectrum from its pdf

How might one generate a sequence of N pseudo random numbers that simultaneously have an e.g. pink spectrum and are gaussian distributed?

A hunch would be to first generate an uncorrelated white uniform pdf sequence, do nonlinear processing to get the desired pdf, then filter it to produce the desired spectrum. But that would alter the pdf. Perhaps the output could be nonlineary processed again, and this could be done iteratively to converge (or not) close to the desired result. But shurely there must be some more solid framework linking the two?

edit:

What if you want an arbitrary PDF and an arbitrary spectrum? Say that I have a random sequence:

N = 1000;
x = 2*rand(N,1) - 1;


What if I want to change this into a pink spectrum, triangular pdf sequence? Or (an approximation of) the spectrum of a Beatles song with a binary pdf

The case for generating a Gaussian distribution is simple as that is exactly what would happen after filtering. A Gaussian distribution of white noise samples is the one case where the spectrum can be modified and the distribution will still be Gaussian. Given any other distribution of white noise samples, once you filter it, the resulting pdf will approach a Gaussian! Given non-white samples, if we filter it over a long enough span to include independent samples in the summation of the filtering, the result will be Gaussian!

As suggested by the central limit theorem, the result after filtering will be Gaussian so it is only the total span which need be adjusted to compensate for the reduced signal power from the filtering. It is for this reason that quantization noise at the output of an FIR filter is Gaussian distributed (which otherwise starts with a uniform distribution). Similarly any white noise distribution will converge to a Gaussian after filtering (as the sum of IID RV's).

As a simple and visual demonstration, consider the case of a uniform distribution after processing with a moving average. First we generate $$2^18$$ uniformly distributed samples and observe the histogram:

test = rand(2**18);
hist(test, 512);


And then we average the samples with a 512 sample moving average filter and observe the histogram after filtering, which is essentially Gaussian (even if this was done over a far shorter average in terms of number of samples, it will still appear quite Gaussian):

>> result = filter(ones(512,1),512,test);
>> hist(result, 512)


As for the more complicated case of going the other way, generating a uniform or any other distribution of a filtered white noise process, please see this answer already provided in this post:

Band-limited random signal with arbitrary distribution?

• Thank you, Dan. That link at the bottom seems very relevant to my (revised) question Apr 16, 2022 at 15:04
• I'm really just picking nits here, but a filtered signal tends toward Gaussian if it's effectively averaged over many samples. So there are edge cases (e.g. $H(z) = 1 - z^{-1}$), where there's not enough samples, and the central limit theorem doesn't apply well. Apr 16, 2022 at 19:20
• @TimWescott yes I think that is a valid point — I originally had written a moving average specifically but then considered that for most practical filtering cases there will be enough delayed samples at about the same level adding resulting in the IID constraint. As soon as you have more than six it seems it is Gaussian enough in most cases Apr 16, 2022 at 19:28

That's actually fairly straight forward.

The easiest method is to generate gaussian distributed white noise and apply a pink filter to it. See https://ccrma.stanford.edu/~jos/sasp/Example_Synthesis_1_F_Noise.html

If you want to be more fancy, you can generate it in the frequency domain. a) create a pink magnitude spectrum, b) add a uniformly distributed random phase, c) do an inverse FFT. See code below.

There are two things to consider:

1. Getting a gaussian distribution is fairly easy. Due to the central limit theorem, most operations you do on a signal will make it more Gaussian. That's why the random phase works
2. "True" pink has infinite energy at DC and very high energy at very low frequencies, which often makes it impractical. In my example I work around this by setting a "hinge point". It's pink down to 20Hz but flat below.

For example

%% script to generate pink noise with a low frequency hinge point
fs = 48000; % sample rate in Hz
fc = 20; % hinge point in Hz
nx = 2^16;  % length of the sequence
ifc = round(fc/fs*nx); % index of the hinge point
% create white magnitude  spectrum
fx = ones(nx/2+1,1);
fx(2:end) = sqrt(ifc./(1:nx/2))';  %  pink roll off excluding DC
fx(1:(ifc+1)) = 1; % set the low end to flat
fx(2:end-1) = fx(2:end-1).*exp(1i*2*pi*rand(nx/2-1,1));
% make symmetric
fx = [fx; conj(fx(end-1:-1:2))];
% inverse FFT and normalize to unity power
pinkNoise = real(ifft(fx));
pinkNoise = pinkNoise./rms(pinkNoise);

% show histograms of pink and white Gaussian noise
clf;
histogram(pinkNoise);
whiteNoise = randn(nx,1);
hold on

histogram(gca,whiteNoise);
legend('Pink','White');
grid('on');


(Dis-)connecting a stochastic signals spectrum from its pdf

There is nothing in particular that connects a signal's spectrum to its PDF. White noise can be Gaussian, it can be binary, it can be pretty much anything you want it to be. Gaussian PDF's don't need to be white, they can have any spectrum you like (within reason)

The main exceptions are sequences that are very short in one domain. Narrow band signals approach the PDF of a sine wave, very short time signals tend to have a whit-ish spectrum.

• Thanks Hilmar. That answers my question very well. Unfortunately, I should have avoided being so specific (see my edited question). Apr 16, 2022 at 15:05

I think that I can answer myself here.

Encouraged by the post provided by Dan Boschen (Band-limited random signal with arbitrary distribution?), I did what I suggested as a "hunch" in my question. It is not generic in that it is not equalizing the histogram, rather going for simpler hard-coded knowledge of the target pdf. But it seems to produce the desired results, and I guess that one could work on a general numerical approximation to the pdf with a little effort.

I still think that this approach is inelegant, which kind of surprise me for such a simple and fundamental task.

%% define target response
%hammock spectrum
K_H = 512;
H = (1+cos(linspace(0, 2*pi, K_H)))/2;
%assume that a binary pdf (-1/+1) of equal probability is required

%% generate white uniform zero-avg pdf input
N = 48*K_H;
x = 2*rand(N,1)-1;

for it = 1:3
%% modify spectrum magnitude
win = sqrt(hann(K));
xb = buffer(x, K, K/2, "nodelay");
X = fft(win.*xb);
P = mean(abs(X), 2);
gs = H' ./ P;
Y = X .* gs;
yb = real(ifft(Y));
x = debuffer(win.*yb, K/2);

%% modify histogram
x(x<0) = -1;
x(x>=0) = 1;
end