# Pink noise and averaging

Disclaimer: I am aware this is somehow at the boundary between physics and statistics, but I have the impression that it is more likely that somebody doing/studying physics, rather than statistics, has stumbled in this problem, which in the end is very practical. Feel free to comment and/or suggest a better community if you disagree!

When measuring any physical parameter, uncertainty can be reduced by repeating the measurement many times and then taking the average value $$\mu=(\sum_1^N x_i)/N$$, or also by measuring for a longer time. Assuming the single measurement has a noise $$\sigma_{x_i}=\sigma_0$$, it is common knowledge that the standard deviation of the average of $$N$$ measurement should scale as $$\sigma_\mu = \sigma_0/\sqrt{N}$$ (or actually $$\sigma_0/\sqrt{N-1}$$)... so the more you measure, the more precise your estimate.

However, indeed this is only true if the $$N$$ measurements are uncorrelated, which is equivalent to stating that the measurement is affected only by white noise. When flicker noise is present (i.e., ultimately, always), the power spectral density of noise deviates from flat in the low-frequency limit, and at some point starts growing as $$\propto 1/\omega$$ (or even steeper but let us neglect that). As a result, this means that the value you measure "now" has some correlation with those you got "previously", and the rule above fails.

Question: How does the standard deviation scale with the number of measurements or with the integration time in the presence of flicker/pink noise?! Can you point me in the right direction and/or do you have a good textbook/reference discussing this? I am having a hard time finding one... statistics books even rarely discuss power spectral densities, not to mention specific cases like this one; in physics texts, it is uncommon that correlated noise is discussed in detail.

What I have concluded so far. It is quite obvious that in the presence of pink noise the precision of the estimate won't improve forever like $$\approx 1/\sqrt{N}$$ by simply increasing the number of measurements $$N$$. I have found vague comments about the existence of a "Flicker floor" but no real clear discussion about it. I see that the problem has a well-defined frequency scale, which is the frequency below which flicker noise takes over and dominates the spectral density, the "Flicker corner" $$\omega_c$$. So I would intuitively expect that measuring for a time that is longer or shorter than the time scale $$1/\omega_c$$ makes a qualitative difference.

One obvious approach to see what happens if I integrate my measurement for a time $$T$$ might be to consider the measurement window - say, $$[-T/2,+T/2]$$ - and use its Fourier transform (sinc function) to see which are the frequencies that I am sampling and calculate the total noise. With white noise, it works fine but when you turn to pink, the integral of the PSD clearly diverges. I suspect this has to do with the fact that even the average does not converge with flicker... but... not sure where to head to from here.

 I have found some interesting clues, "flicker floor" was the main search key that led me to the link below. I have not gone through the work in details but it seems to answer to at least some of my doubts

https://arxiv.org/abs/1407.7760v6

• Very short comment: if you apply the PSD of the pink noise to use when "generating" samples, you'll start to see how the number of samples needs to be "weighted" related to the probability curve. Commented Jun 5, 2021 at 14:05
• You may have to analyze a covariance matrix of all measurements, rather than their standard devuations (diagonal elements of this matrix). Also, it might be worthwhile to look at it as time series (using one of the available models: ARMA, ARIMA, etc.) Just to give you some keywords. Commented Jun 5, 2021 at 14:27
• Thanx Roger and Carl! I try and look further into the thing with these inputs. Indeed I tried once to compute the covariance matrix elements and proceed from there: I cant remeber how it ended exactly, I just recall I got lost but I might try again.
– Ste
Commented Jun 5, 2021 at 16:49
• @Ste there's a basic problem that the noise power is infinite if the 1/f power spectral density extends all the way down to zero. Nevertheless, if the wanted signal has a distinct spectrum (i.e. not simply a dc constant value) , there are methods to optimally estimate (filter) the signal. Commented Jun 11, 2021 at 4:00

Secondly, without using a simulation, as you mentioned, for independent sequence of variables, of strict sence stationary process~(SSS), there is a Central Limit Theorem (CLT), which says that the distribution $$X_n = \sum_{k=1}^{n} \frac{x[k]}{n}$$, and it converges to $$\bar{X}_n \sim N(0,\sigma \sqrt{n})$$~(by the way I am going to consider only centered process, it is easy to generalize for non-centered one). It means that taking n samples of the process the average is an expected value of: $$E\left\{ \sum_{k=1}^{n} \frac{x[k]}{n} \right\} = \sum_{k=1}^{n} \frac{0}{n} = 0$$ Simultaneously, the standard deviation of this distribution is $$\sigma^2_{\bar{X}_n} = E\left\{ \bar{X}^2_n \right\} = E\left\{ \left( \sum_{k=1}^{n} \frac{x[k]}{n} \right)^2\right\}$$ So it is necessary to square the sum and the only surviving terms due to the centered propriety are $$E\{ x[k] x[k] \} = \sigma^2$$ $$\sigma^2_{\bar{X}_n} = \frac{n \sigma^2}{n^2} = \frac{\sigma^2}{n}$$
Next stage is to say that process you describe is SSS, has a certain distribution and it is pink. The last means that, as you mentioned, that there is a correlation between closely spaced time samples. So, first of all the average value of the $$E\{\bar{X}_n\} = 0$$, since expected value is a linear operation and each sample has zero expected value. The standard deviation in this case can be estimated again: $$\sigma^2_{\bar{X}_n} = E\left\{ \bar{X}^2_n \right\} = E\left\{ \left( \sum_{k=1}^{n} \frac{x[k]}{n} \right)^2\right\}$$ So, once again it is necessary to square all the terms, where we get: $$\sigma^2_{\bar{X}_n} = \frac{n \sigma^2 + 2 E\{ x[1]x[2]\} + 2 E\{ x[1]x[3] \} + ... }{n^2}$$ So covariances of different samples do appear. Now, since a process SSS, its auto correlation is $$R[n,m] = E\{ x[n] x[m] \} = R[m-n]$$ Thus, we keep only terms: $$\sigma^2_{\bar{X}_n} = \frac{n \sigma^2 + 2 R[1] (n-1) + 2 R[2] (n-2) + ... + R[n-1] }{n^2} = \frac{\sigma}{n} + 2R[1] \frac{n-1}{n^2} + 2R[2] \frac{n-2}{n^2} +... + 2R[n-1]\frac{1}{n^2}$$ So there is different weight of each autocorrelation term, and for example, les weight has $$R[n-1]$$ since, evidently it appears only once. Also, it is important, that if there is no covariance between sample (= white noise), we get $$\sigma^2_{\bar{X}_n} = \frac{\sigma^2}{n}$$ as in CLT.