General system--windowed estimation of correlation spectrum

I assume that ${\mathcal{F}}(t)$ (which I call the `signal') for a single, ergodic trajectory has been generated and stored at a fine enough time-resolution $t_r$ to capture all the desired features. This can be done by evaluating $-\partial {\mathcal{H}}/\partial x$ along the trajectory, and subtracting the average value (estimated, or often analytically known). We need $t_r < \omega_0$, where $\omega_0$ is the upper frequency limit desired for the band profile.

The auto-correlation spectrum of this single signal is $\tilde{C}_{\mathcal{F}}(\omega)$, and if the signal exists over all times $-\infty < t < \infty$ then $\tilde{C}_{\mathcal{F}}(\omega)$ is a stochastic quantity uncorrelated in $\omega$. The microcanonical average $\tilde{C}(\omega)= \langle \tilde{C}_{\mathcal{F}}(\omega)\rangle_{{\mbox{\tiny E}}}$ (I drop the energy subscript) is reached by averaging over initial conditions of trajectories. It is $\tilde{C}(\omega)$ that is the desired `band profile'. However, $\tilde{C}(\omega)$ is also the local mean of $\tilde{C}_{\mathcal{F}}(\omega)$, so can be more easily estimated by smoothing $\tilde{C}_{\mathcal{F}}(\omega)$. This smoothing (convolution with Gaussian of width $\omega_s$) is performed using the Fast Fourier Transform (FFT) [161]. Since computationally the signal must be of necessarily finite length $\sim t_0$, correlations in $\omega$-space arise in $\tilde{C}_{\mathcal{F}}(\omega)$ on the scale $t_0^{-1}$. Therefore the smoothing operation corresponds to averaging $N \approx \omega_s t_0$ independent samples, and a fractional error of $N^{-1/2}$ results for the band profile estimate. This estimation error appears as multiplicative noise with frequency correlation scale $\omega_s$. Typical parameters were $\omega_0 = 4, \omega_s = 10^{-2}$, $t_0 = 10^5$, giving about 3% RMS estimation error.

There is another `failure mode' of this estimation procedure: There are only $M = \sim t_0/\tau_{{\mbox{\tiny cl}}}$ independent samples in the finite sample length, therefore only $M$ independent values can appear in $cfw$. If $M \ll \omega_0 t_0$ then the errors may exceed that quoted above, because samples are no longer independent. This failure mode only arises when $\tilde{C}(\omega)$ is desired at large $\omega$, which forces $t_r$ to shrink, and therefore $t_0$ to also shrink if the memory requirements are not to grow. In this case, true averaging over trajectories is the only remedy.

Now I describe how to find $\tilde{C}_{\mathcal{F}}(\omega)$ given ${\mathcal{F}}(t)$. I will always use the Fourier Transform (FT) convention

$$\tilde{A}(\omega) \; = \; \int_{-\infty}^{\infty} \!\!dt \, e^{i \omega t} A(t) \hspace{1.4in} \mbox{forward},$$

$$A(t) \; = \; \frac{1}{2\pi} \int_{-\infty}^{\infty} \!\!d\omega \, e^{-i \omega t} \tilde{A}(\omega) \hspace{1in} \mbox{backward}.$$ (B.1)

For this general real signal $A(t)$, the Wiener-Khinchin Theorem (easy to prove [161]) equates the power spectrum to the FT of the auto-correlation,

$$\vert\tilde{A}(\omega) \vert^2 \;= \; \tilde{C}_A(\omega) ,... ...; \equiv \; \int_{-\infty}^{\infty} \!\!dt \, A(t) A(t+\tau).$$ (B.2)

I choose $A(t)$ to be a windowed version of the signal,

$$A(t) \; = \; W(t) \, {\mathcal{F}}(t), \hspace{1in} \mbox{$W(t) = $\ windowing function},$$ (B.3)

where $W(\pm\infty)\rightarrow0$ and has finite characteristic width $t_0$. This windowing function defines the maximum trajectory length needed. Therefore the auto-correlation is related to the desired $C_{\mathcal{F}}(\tau)$ by

$$C_A(\tau) \; = \; \int_{-\infty}^{\infty} \!\!dt \, W(t) W... ...\mathcal{F}}(t+\tau) \; = \; C_W(\tau) C_{\mathcal{F}}(\tau),$$ (B.4)

where the last step involved identifying $W(t) W(t+\tau)$ as a $\tau$-dependent weighting function whose $t$-integral is $C_W(\tau)$, giving a generalization of the `top-hat' weighting function of e.g. Eq.(2.30). Applying the Convolution Theorem[161] and dividing by $Z$,

$$\frac{\vert\tilde{A}(\omega)\vert^2}{Z} \; = \; \frac{1}{Z}(\tilde{C}_{\mathcal{F}} \ast \tilde{C}_W)(\omega) .$$ (B.5)

$\tilde{C}_W(\omega)$ is a narrowly peaked function $\approx Z \delta(\omega)$ where the normalization is $Z = 2\pi \int_{-\infty}^{\infty} \!\!dt \, W(t)^2$. Therefore we have the LHS as our estimate for $\tilde{C}_{\mathcal{F}}(\omega)$. Structure below frequency scale $t_0^{-1}$ (the $\delta$-function width) is lost. I typically used the Gaussian window $W(t) = e^{-t^2/(2t_0^2)}$, for which $Z = 2 \pi^{3/2} t_0$. The maximum signal length $T_0$ needed to evaluate the windowed signal $A(t)$ was typically chosen as $10 t_0$, beyond which $W(t)$ is negligible. Because $W(t)$ dies to zero very smoothly, the windowed signal $A(t)$ of length $T_0$ can now be used as a periodic discretely-sampled array, and the FFT used to find $\tilde{A}(\omega)$.

All FFTs were implemented on a Compaq XP1000 (21264 667MHz Alpha processor) running C++ and the Compaq Extended Math Library (CXML). Using this setup, an 8 million ($2^{23}$) point double-precision FFT takes about 5 seconds.