Figure 1

A realization of the random process $x\left(t\right)$ [Eq. (2)] with $A=1$, $n=1$, and with fixed decay rate $\lambda =0.5$ (all quantities are given in arbitrary units; $\lambda$ is given in inverse time units). The single exponential decays are clearly visible, and the random process is obviously non-Gaussian.Reuse & Permissions

Figure 2

Structure and update dynamics of the linked list that holds the Poisson-distributed transition events. (a) Structure of the list at the $j\text{th}$ sampling time $s_j$; each node contains a variable that points to the next node (the end of the list is marked by a null pointer), and stores the time $t_k$ of the transition and the decay rate ${\lambda }_k$ of the associated pulse response function. The list contains only nodes such that $s_j-t_k⩽N_{\mathit{\text{decay}}}∕{\lambda }_k$. The response of the system is computed from the sum ${\sum }_k\exp [-{\lambda }_k(s_j-t_k)]$, where the index $k$ ranges over all the list elements such that $s_j-t_k⩾0$ (the list head is usually excluded). (b) If the $(j+1)\text{th}$ sampling time is greater than the time stored in the list head (as is usually the case), the program generates as many transition times as needed to reach (and possibly overcome) the $(j+1)\text{th}$ sampling time (light gray boxes in the figure, primed quantities), and next it scans the list to discard all the nodes such that $s_{j+1}-t_k>N_{\mathit{\text{decay}}}∕{\lambda }_k$ (dark gray boxes). At this point the program computes the new response and steps to the next sampling step.Reuse & Permissions

Figure 3

Length of the linked list in a simulation with $A=1$ and a single decay rate $\lambda =0.001$: the linked list is initially empty, and it fills up at a constant rate. In this case $n=1$, $N_{\mathit{\text{decay}}}=20$, and $\Delta s=1$ and therefore the fill-up time is $(N_{\mathit{\text{decay}}}∕\lambda )=20000$ the fill-up length is $n(N_{\mathit{\text{decay}}}∕\lambda )=20000$, and the number of samples required for the initial fill-up is $(N_{\mathit{\text{decay}}}∕\lambda )∕\Delta s=20000$. After the initial fill-up the length of the linked list fluctuates about the average filling level.Reuse & Permissions

Figure 4

Plot of the normalized signal amplitude $[x\left(t\right)-⟨x⟩]∕\sigma$ [$\sigma$ is the standard deviation of the amplitude, i.e., the square root of the variance (6)] in the simulation run described in Fig. 3 and in the text. At the beginning the linked list which contains the process memory is empty, and the signal is very far off the predicted average; as the list fills up to level, the signal quickly reaches the predicted average.Reuse & Permissions

Figure 5

Detail of the normalized signal amplitude shown in Fig. 4, just after the list has filled up to the average level. This signal displays the large characteristic upward an downward swings that are well known in the theory of random walks [23].Reuse & Permissions

Figure 6

The histogram shows the amplitude distribution of 262 144 samples from the realization of the random process $x\left(t\right)$ shown in Fig. 4, after the list fill-up. The continuous curve is a Gaussian with the mean and standard deviation estimated from the samples. There is no visible skewness, because in this simulation run $\lambda ∕n=0.001$, which corresponds to a very low skewness (23), but there are multiple peaks, which are due to the nonstationarity of a true $1∕f^2$ process (which is well approximated here), and which require an extremely long observation time to establish the Gaussianity of the process [23].Reuse & Permissions

Figure 7

DFT spectrum obtained from 262 144 samples from the realization of the normalized signal amplitude $[x\left(t\right)-⟨x⟩]∕\sigma$ shown in Fig. 4. The continuous curve shows the theoretical power spectral density (10). Because of sampling without low-pass filtering there is some aliasing and the DFT spectrum shows a slight upward bend at high frequency. Since the sampling interval is $\Delta s=1$, the Nyquist (angular) frequency is just ${\omega }_{\mathit{\text{Nyquist}}}=\pi$ (here and in the following spectra time is measured in arbitrary units as in the previous figures, and frequency units are defined accordingly). The arrow marks the position of the single decay rate in this simulation $\lambda =0.001$.Reuse & Permissions

Figure 8

Averaged DFT spectrum obtained from the same 262 144 samples as the spectrum in Fig. 7, split into 32 blocks of 8192 samples each. The continuous curve shows the theoretical power spectral density (10), and now it includes also the first-order correction to aliasing. Because of the low-frequency correlation between the blocks (that have been obtained from the same simulation record), the average spectrum is a bit higher than expected and the theoretical prediction has been globally shifted 20% higher to fit the average spectrum; this artifact is absent in the analysis of the whole record (the low-frequency plateau of the spectrum in Fig. 7 fits the theoretical curve exactly as expected). As in Fig. 7, the arrow marks the position of the single decay rate in this simulation $\lambda =0.001$: because of the shorter record length used for DFT analysis, the frequency resolution is poorer here, and the spectrum mimics quite well the behavior of a true $1∕f^2$ spectrum, over about three frequency decades.Reuse & Permissions

Figure 9

Length of the linked list in a simulation with $A=1$ and a uniform distribution of decay rates in the range ${\lambda }_{\mathit{min}}=0.0001$, ${\lambda }_{\mathit{max}}=1$: the linked list is initially empty, at it fills up with a variable rate that depends on the distribution of decay rates. In this case $n=10$, $N_{\mathit{\text{decay}}}=20$, $\Delta s=1$, and $⟨1∕\lambda ⟩\approx 9.211$, and therefore the fill-up time is $N_{\mathit{\text{decay}}}∕{\lambda }_{\mathit{min}}=200000$ the fill-up length is $n\left(N_{\mathit{\text{decay}}}⟨1∕\lambda ⟩\right)\approx 1842$, and the number of samples required for the initial fill-up is $(N_{\mathit{\text{decay}}}∕{\lambda }_{\mathit{min}})∕\Delta s=200000$. After the initial fill-up the length of the linked list fluctuates about the average filling level.Reuse & Permissions

Figure 10

Detail of the normalized signal amplitude in the simulation of Fig. 9, just after the list has filled up to the average level.Reuse & Permissions

Figure 11

The histogram shows the amplitude distribution of 262 144 samples from the realization of the random process $x\left(t\right)$ shown in Fig. 10, after the list fill-up. The continuous curve is a Gaussian with the mean and standard deviation estimated from the samples. In contrast to the histogram in Fig. 6, now the amplitude distribution appears slightly skewed, because in this simulation run $1∕\left(n⟨1∕\lambda ⟩\right)\approx 0.0109$, noticeably higher than the corresponding value for Fig. 6.Reuse & Permissions

Figure 12

Averaged DFT spectrum obtained from the 262 144 samples in the simulation of Fig. 9, split into 32 blocks of 8192 samples each. The continuous curve shows the theoretical power spectral density (13), which includes also the first-order correction to aliasing. The arrows mark the positions of the extreme decay rates ${\lambda }_{\mathit{min}}=0.0001$ and ${\lambda }_{\mathit{max}}=1$. The spectral resolution $\Delta \omega \approx 0.0015$ is larger than the minimum decay rate ${\lambda }_{\mathit{min}}$, and the spectrum mimics quite well the behavior of a true $1∕f$ spectrum (dashed line), over more than three frequency decades.Reuse & Permissions

Figure 13

Averaged DFT spectrum obtained from $2^{20}=1048576$ samples with $A=1$ and a range-limited power-law distribution of decay rates ${\lambda }^{-\beta }$, with $\beta =0.2$, in the range ${\lambda }_{\mathit{min}}=0.0001$, ${\lambda }_{\mathit{max}}=1$, split into 32 blocks of 32 768 samples each. The continuous curve shows the theoretical power spectral density (16), which includes also the first-order correction to aliasing. The arrows mark the positions of the extreme decay rates ${\lambda }_{\mathit{min}}=0.0001$ and ${\lambda }_{\mathit{max}}=1$. The spectral resolution $\Delta \omega \approx 0.00038$ is larger than the minimum decay rate ${\lambda }_{\mathit{min}}$, and the spectrum mimics quite well the behavior of a true $1∕f^{1.2}$ spectrum (dashed line), over more than three frequency decades.Reuse & Permissions

Figure 14

Ratio $\mathrm{var}\left[\mathrm{Re}\left(F_k\right)\right]∕S_k$ obtained from $2^{20}=1048576$ samples with $A=1$ and a uniform distribution of decay rates, in the range ${\lambda }_{\mathit{min}}=0.0001$, ${\lambda }_{\mathit{max}}=1$, averaged over the DFT results obtained from 128 blocks of 8192 samples each. The average ratio fluctuates about a constant level, and this is in line with the usual hypothesis that $\mathrm{var}\left[\mathrm{Re}\left(F_k\right)\right]\propto S_k$ (see, e.g., [8]).Reuse & Permissions