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

The auto-correlation spectrum of this single signal is
, and if
the signal exists over all times
then
is a stochastic quantity uncorrelated in .
The microcanonical average
(I drop the energy subscript)
is reached by averaging over initial conditions of
trajectories.
It is
that is the desired `band profile'.
However,
is also the *local mean* of
, so can be more easily
estimated by smoothing
.
This smoothing (convolution with Gaussian of width )
is performed using the Fast Fourier Transform (FFT) [161].
Since computationally the signal must be of necessarily finite length
,
correlations in -space arise in
on the scale .
Therefore the smoothing operation corresponds to averaging
independent samples, and a fractional error of
results for the band profile estimate.
This estimation error appears as multiplicative noise with
frequency correlation scale .
Typical parameters were
, , giving
about 3% RMS estimation error.

There is another `failure mode' of this estimation procedure: There are only independent samples in the finite sample length, therefore only independent values can appear in . If then the errors may exceed that quoted above, because samples are no longer independent. This failure mode only arises when is desired at large , which forces to shrink, and therefore 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
given
.
I will always use the Fourier Transform (FT) convention

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

I choose to be a windowed version of the signal,

where and has finite characteristic width . This windowing function defines the maximum trajectory length needed. Therefore the auto-correlation is related to the desired by

where the last step involved identifying as a -dependent weighting function whose -integral is , giving a generalization of the `top-hat' weighting function of

is a narrowly peaked function where the normalization is . Therefore we have the LHS as our estimate for . Structure below frequency scale (the -function width) is lost. I typically used the Gaussian window , for which . The maximum signal length needed to evaluate the windowed signal was typically chosen as , beyond which is negligible. Because dies to zero very smoothly, the windowed signal of length can now be used as a periodic discretely-sampled array, and the FFT used to find .

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 () point double-precision FFT takes about 5 seconds.