StocProc_FFT

class stocproc.stocproc.StocProc_FFT(spectral_density, t_max, bcf_ref, intgr_tol=0.01, intpl_tol=0.01, seed=None, negative_frequencies=False, scale=1)

Simulate Stochastic Process using FFT method

This method uses the relation of the auto correlation to the non negative real valued spectral density \(J(\omega)\). The integral can be approximated by a discrete integration scheme

\[\alpha(\tau) = \int_{\omega_\mathrm{min}}^{\omega_\mathrm{max}} \mathrm{d}\omega \, \frac{J(\omega)}{\pi} e^{-\mathrm{i}\omega \tau} \approx \sum_{k=0}^{N-1} w_k \frac{J(\omega_k)}{\pi} e^{-\mathrm{i} k \omega_k \tau}\]

where the weights \(\omega_k\) depend on the particular integration scheme. For a process defined as

\[Z(t) = \sum_{k=0}^{N-1} \sqrt{\frac{w_k J(\omega_k)}{\pi}} Y_k \exp^{-\mathrm{i}\omega_k t}\]

with independent complex random variables \(Y_k\) such that \(\langle Y_k \rangle = 0\), \(\langle Y_k Y_{k'}\rangle = 0\) and \(\langle Y_k Y^\ast_{k'}\rangle = \Delta \omega \delta_{k,k'}\) it is easy to see that its auto correlation function will be exactly the approximated auto correlation function.

\[\begin{split}\begin{align} \langle Z(t) Z^\ast(s) \rangle = & \sum_{k,k'} \frac{1}{\pi} \sqrt{w_k w_{k'} J(\omega_k)J(\omega_{k'})} \langle Y_k Y_{k'}\rangle \exp(-\mathrm{i}(\omega_k t - \omega_k' s)) \\ = & \sum_{k} \frac{w_k}{\pi} J(\omega_k) e^{-\mathrm{i}\omega_k (t-s)} \approx & \alpha(t-s) \end{align}\end{split}\]

To calculate \(Z(t)\) the Discrete Fourier Transform (DFT) can be applied as follows:

\[Z(t_l) = e^{-\mathrm{i}\omega_\mathrm{min} t_l} \sum_{k=0}^{N-1} \sqrt{\frac{w_k J(\omega_k)}{\pi}} Y_k e^{-\mathrm{i} 2 \pi \frac{k l}{N} \frac{\Delta \omega \Delta t}{ 2 \pi} N}\]

Here \(\omega_k\) has to take the form \(\omega_k = \omega_\mathrm{min} + k \Delta \omega\) and \(\Delta \omega = (\omega_\mathrm{max} - \omega_\mathrm{min}) / (N-1)\) which limits the itegration schemes to those with equidistant weights. For the DFT scheme to be applicable \(\Delta t\) has to be chosen such that \(2\pi = N \Delta \omega \Delta t\) holds. Since \(J(\omega)\) is real it follows that \(X(t_l) = X^\ast(t_{N-l})\). For that reason the stochastic process has only \((N+1)/2\) (odd \(N\)) and \((N/2 + 1)\) (even \(N\)) independent time grid points.

To generate a process with given auto correlation function on the interval [0, t_max] requires that the auto correlation function approximation is valid for all t in [0, t_max].

This is ensured by automatically determining the number of sumands N and the integral boundaries \(\omega_\mathrm{min}\) and \(\omega_\mathrm{max}\) such that discrete Fourier transform of the spectral density matches the desired auto correlation function within the tolerance intgr_tol for all discrete \(t_l \in [0, t_\mathrm{max}]\).

As the time continuous process is generated via cubic spline interpolation, the deviation due to the interpolation is controlled by the parameter intpl_tol. The maximum time step \(\Delta t\) is chosen such that the interpolated valued at each half step \(t_i + \Delta t /2\) differs at most intpl_tol from the exact value of the auto correlation function.

If not fulfilled already N and the integration boundaries are increased such that the \(\Delta t\) criterion from the interpolation is met.

Parameters:

• spectral_density – the spectral density \(J(\omega)\) as callable function object
• t_max – \([0,t_\mathrm{max}]\) is the interval for which the process will be calculated
• bcf_ref – a callable which evaluates the Fourier integral exactly
• intgr_tol – tolerance for the integral approximation
• intpl_tol – tolerance for the interpolation
• seed – if not None, use this seed to seed the random number generator
• negative_frequencies – if False, keep \(\omega_\mathrm{min} = 0\) otherwise find a negative \(\omega_\mathrm{min}\) appropriately just like :math:`omega_mathrm{max}

Todo

implement bcf_ref = None and use numeric integration as default

calc_z(y)

calculate

\[Z(t_l) = e^{-\mathrm{i}\omega_\mathrm{min} t_l} \mathrm{DFT}\left( \sqrt{\frac{w_k J(\omega_k)}{\pi}} Y_k \right)\]

and return values with \(t_l < t_\mathrm{max}\)

static get_key(t_max, bcf_ref, intgr_tol=0.01, intpl_tol=0.01)

get_num_y()

The number of independent random variables Y is given by the number of discrete times \(t_l < t_\mathrm{max}\) from the Fourier Transform

get_time()

Returns the time \(t_k\) corresponding to the values \(z_k\)

These times are determined by the integration weights method.

get_z()

Returns the discrete process \(z_k\).

new_process(y=None, seed=None)

generate a new process by evaluating calc_z() with new random variables \(Y_i\)

Parameters:

• y – independent normal distributed complex valued random variables with \(\sigma_{ij}^2 = \langle y_i y_j^\ast \rangle = 2 \delta_{ij}\)
• seed – if not None set seed to seed before generating samples

When y is given use these random numbers as input for calc_z() otherwise generate a new set of random numbers.

set_scale(scale)