Obtaining 1/f noise from superposition of Lorentzians
Earlier we have taken a look at a power spectral density of a signal generated by a single charge carrier. We have seen that it generates a signal composed of pulses and gaps. If the pulse and gap duration is distributed according to the exponential distribution (i.e., Poisson process model applies), then the power spectral density of the signal has a characteristic Lorentzian shape. Which is nothing alike 1/f noise, which we are looking for. But maybe we can still obtain it?
Superposition of Lorentzians
So, earlier we have obtained that the power spectral density generated by a single charge carrier, for which the Poisson process model applies, is given by
\begin{equation} S\left(f\right)=\frac{4a^{2}\bar{\nu}}{\left(\lambda_{\theta}+\lambda_{\tau}\right)^{2}+4\pi^{2}f^{2}}. \end{equation}
Here for the sake of simplicity let us assume that \( \lambda_{\theta} = \lambda_{\tau} \). As the rates are assumed to be equal let us drop the indices. This allows us to a simpler form for the power spectral density in our further derivations, i.e.,
\begin{equation} S\left(f\right)=\frac{a^{2}\bar{\nu}}{\lambda^{2}+\pi^{2}f^{2}}. \end{equation}
Our model will rely on superposition of these Lorentzian shape with different characteristic rates \( \lambda \). As is well known, the Fourier transform of a sum of two signals (functions) is the same as the sum of their transforms:
\begin{equation} \mathcal{F}\left\{f(t)+g(t)\right\} =\mathcal{F}\left\{f(t)\right\}+\mathcal{F}\left\{g(t)\right\}. \end{equation}
Power spectral density of the sum of two signals (functions) is not necessarily a sum of their individual power spectral densities. If the two signals are correlated, there will be cross spectral density term. But if the two signals are independent (as we assume is the case here), we will have just an ordinary sum of the components.
Let us assume that the distribution of characteristic rates is continuous uniform distribution (with certain \( \lambda_{\text{min}} \) and \( \lambda_{\text{max}} \) bounds). Then:
\begin{equation} S\left(f\right)=\frac{1}{\lambda_{\text{max}}-\lambda_{\text{min}}}\int_{\lambda_{\text{min}}}^{\lambda_{\text{max}}} \frac{a^{2}\bar{\nu}}{\lambda^{2}+\pi^{2}f^{2}} d \lambda = \frac{a^2 \bar{\nu}\left[\operatorname{arccot}\left(\frac{\pi f}{\lambda_{\text{max}}}\right)-\operatorname{arccot}\left(\frac{\pi f}{\lambda_{\text{min}}}\right)\right]}{(\lambda_{\text{max}}-\lambda_{\text{min}}) \pi f}. \end{equation}
At first there doesn't seem to be any problem with this, but recall that \( \bar{\nu} \) represents mean number of pulses per unit time. So, it does depend on \( \lambda \). In our particular case (when the pulse and gap duration distributions match) we have that \( \bar{\nu} = \lambda / 2 \). To quickly fix this let us dynamically scale the pulse height
\begin{equation} a = \sqrt{\frac{1}{\bar{\nu}}} = \sqrt{\frac{2}{\lambda}}. \end{equation}
Then we simply have that \( a^2 \bar{\nu} = 1 \) for any \( \lambda \), and we can ignore this term as we already did in the derivation above.
Interactive app
To effectively show the key point of this post (i.e., that you can get 1/f noise out of superposition of Lorentzian shapes) here we do not generate the signals themselves. Unlike in the earlier posts here instead we just add the Lorentzian shapes, which should emerge as a result of long duration detailed simulations. Simulating individual signals would likely just take unreasonably long until properly looking result would emerge.
So, in the app below you can adjust the bounds of the \( \lambda \) distribution. In the plot the simulation results are shown as red dots, and they quickly converge towards analytically derived formula (shown as black curve). As you change the bound of the \( \lambda \) distribution observe how the range of frequencies over which 1/f noise is observed changes. Note that for the lowest frequencies white noise dominates, while for the highest frequencies \( S(f) \sim 1/f^2 \) dependence emerges.