diff --git a/docs/source/StocProc_ABC.rst b/docs/source/StocProc_ABC.rst
new file mode 100644
index 0000000..30c3127
--- /dev/null
+++ b/docs/source/StocProc_ABC.rst
@@ -0,0 +1,8 @@
+Abstract Stochastic Process Class
+---------------------------------
+
+.. py:currentmodule:: stocproc.stocproc
+.. autoclass:: stocproc.stocproc._abcStocProc
+    :members:
+
+    .. automethod:: __call__
diff --git a/docs/source/StocProc_FFT.rst b/docs/source/StocProc_FFT.rst
index fa16ad9..433ec69 100644
--- a/docs/source/StocProc_FFT.rst
+++ b/docs/source/StocProc_FFT.rst
@@ -1,12 +1,9 @@
 StocProc_FFT
 ------------
 
+.. py:currentmodule:: stocproc.stocproc
 .. autoclass:: stocproc.stocproc.StocProc_FFT
     :members:
     :inherited-members:
     :undoc-members:
 
-.. autoclass:: stocproc.stocproc.StocProc_TanhSinh
-    :members:
-    :inherited-members:
-    :undoc-members:
\ No newline at end of file
diff --git a/docs/source/StocProc_KLE.rst b/docs/source/StocProc_KLE.rst
index 2500c38..77c0945 100644
--- a/docs/source/StocProc_KLE.rst
+++ b/docs/source/StocProc_KLE.rst
@@ -10,12 +10,13 @@ StocProc_KLE
 
     .. automethod:: __call__
 
-.. py:currentmodule:: stocproc.method_kle
-.. autofunction:: stocproc.method_kle.solve_hom_fredholm
-.. autofunction:: stocproc.method_kle.get_mid_point_weights_times
-.. autofunction:: stocproc.method_kle.get_trapezoidal_weights_times
-.. autofunction:: stocproc.method_kle.get_simpson_weights_times
-.. autofunction:: stocproc.method_kle.get_four_point_weights_times
-.. autofunction:: stocproc.method_kle.get_gauss_legendre_weights_times
-.. autofunction:: stocproc.method_kle.get_tanh_sinh_weights_times
-.. autofunction:: stocproc.method_kle.auto_ng
+..
+    .. py:currentmodule:: stocproc.method_kle
+    .. autofunction:: stocproc.method_kle.solve_hom_fredholm
+    .. autofunction:: stocproc.method_kle.get_mid_point_weights_times
+    .. autofunction:: stocproc.method_kle.get_trapezoidal_weights_times
+    .. autofunction:: stocproc.method_kle.get_simpson_weights_times
+    .. autofunction:: stocproc.method_kle.get_four_point_weights_times
+    .. autofunction:: stocproc.method_kle.get_gauss_legendre_weights_times
+    .. autofunction:: stocproc.method_kle.get_tanh_sinh_weights_times
+    .. autofunction:: stocproc.method_kle.auto_ng.
diff --git a/docs/source/StocProc_TS.rst b/docs/source/StocProc_TS.rst
new file mode 100644
index 0000000..1902ced
--- /dev/null
+++ b/docs/source/StocProc_TS.rst
@@ -0,0 +1,7 @@
+StocProc_TanhSinh
+-----------------
+
+.. autoclass:: stocproc.stocproc.StocProc_TanhSinh
+    :members:
+    :inherited-members:
+    :undoc-members:
\ No newline at end of file
diff --git a/docs/source/api.rst b/docs/source/api.rst
new file mode 100644
index 0000000..dd7242a
--- /dev/null
+++ b/docs/source/api.rst
@@ -0,0 +1,34 @@
+API documentation
+=================
+
+The ``_abcStocProc`` interface definition
+
+.. toctree::
+   StocProc_ABC
+
+and the classes implementing the ``_abcStocProc`` interface.
+
+.. toctree::
+   :maxdepth: 3
+
+   StocProc_FFT
+   StocProc_KLE
+   StocProc_TS
+
+Function on the mudule level
+----------------------------
+
+.. autofunction:: stocproc.loggin_setup
+.. autofunction:: stocproc.version
+.. autofunction:: stocproc.version_full
+
+Method related submodules
+-------------------------
+
+.. toctree::
+   :maxdepth: 3
+
+   method_kle
+   method_ft
+
+
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 25252e7..119e113 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -68,9 +68,12 @@ author = 'Richard Hartmann'
 # built documents.
 #
 # The short X.Y version.
-version = '1.0'
+
+import stocproc
+
+version = stocproc.version()
 # The full version, including alpha/beta/rc tags.
-release = '1.0.0'
+release = stocproc.version_full()
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
diff --git a/docs/source/example.rst b/docs/source/example.rst
index 59fcf1b..d2d368a 100644
--- a/docs/source/example.rst
+++ b/docs/source/example.rst
@@ -1,5 +1,5 @@
-Example
-=======
+Full Example
+============
 
 This example sets up a generator for stochastic processes with exponential
 auto correlation. A single process is shown as well as the reconstruction
diff --git a/docs/source/index.rst b/docs/source/index.rst
index ff6f87a..6429134 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -1,2 +1,78 @@
-.. automodule:: stocproc
+StocProc
+========
 
+The StocProc module is a Python3 module allowing to sample
+Gaussian stochastic processes which are wide-sense stationary,
+continuous in time and complex valued. In particular for a given auto correlation function
+:math:`\alpha(\tau)` the stochastic process :math:`z(t)` obeys.
+
+.. math:: \langle z(t) \rangle = 0 \qquad \langle z(t)z(s) \rangle = 0 \qquad \langle z(t)z^\ast(s) \rangle = \alpha(t-s)
+
+Here :math:`\langle \cdot \rangle` denotes the ensemble average.
+
+
+Example
+-------
+
+The example will setup a process generator for an exponential auto correlation function
+and sample a single realization.
+
+::
+
+    def lsd(w):
+        # Lorenzian spectral density
+        return 1/(1 + (w - _WC_)**2)
+
+
+    def exp_ac(t):
+        # exponential auto correlation function
+        # note there is a factor of one over pi in the definition of the auto correlation function
+        # exp_ac(t) = 1/pi int_{-infty}^infty d w  lsd(w) exp(-i w t)
+        return np.exp(- np.abs(t) - 1j*_WC_*t)
+
+    _WC_ = 5
+    t_max = 10
+    print("setup process generator")
+    stp = sp.StocProc_FFT(spectral_density = lsd,
+                          t_max = t_max,
+                          bcf_ref = exp_ac,
+                          intgr_tol=1e-2,
+                          intpl_tol=1e-2):
+
+    print("generate single process")
+    stp.new_process()
+    zt = stp()    # get discrete process
+
+.. image:: ../../examples/proc.*
+    :width: 500px
+    :align: center
+    :height: 300px
+
+The full example code can be found :doc:`here </example>`.
+
+
+How to Install
+--------------
+
+Install the latest version via pip ::
+
+    pip install stocproc
+
+or fetch the bleeding edge version from the `git repository <https://github.com/cimatosa/stocproc>`_ ::
+
+    git clone https://github.com/cimatosa/stocproc.git
+
+and install the package invoking ``setup.py``.
+
+Note: The stocproc module depends on `fcSpline <https://github.com/cimatosa/fcSpline>`_ , a fast cubic spline interpolator for equally distributed nodes.
+
+
+Documentation
+-------------
+
+.. toctree::
+   :maxdepth: 3
+
+   methods
+   api
+   example
diff --git a/docs/source/method_ft.rst b/docs/source/method_ft.rst
new file mode 100644
index 0000000..dab60ec
--- /dev/null
+++ b/docs/source/method_ft.rst
@@ -0,0 +1,10 @@
+method_ft Module
+-----------------
+
+helper functions for the FT method
+
+.. py:currentmodule:: stocproc
+
+.. automodule:: stocproc.method_ft
+    :members:
+    :undoc-members:
diff --git a/docs/source/method_kle.rst b/docs/source/method_kle.rst
new file mode 100644
index 0000000..6196dbd
--- /dev/null
+++ b/docs/source/method_kle.rst
@@ -0,0 +1,9 @@
+method_kle Module
+-----------------
+
+helper functions for the KLE method
+
+
+.. automodule:: stocproc.method_kle
+    :members:
+    :undoc-members:
diff --git a/docs/source/methods.rst b/docs/source/methods.rst
new file mode 100644
index 0000000..6a960f4
--- /dev/null
+++ b/docs/source/methods.rst
@@ -0,0 +1,82 @@
+Methods to Sample Stochastic Processes
+======================================
+
+
+This module allows to sample Gaussian stochastic processes which are wide-sense stationary,
+continuous in time and complex valued. In particular for a given auto correlation function
+:math:`\alpha(\tau)` the stochastic process :math:`z(t)` obeys.
+
+.. math:: \langle z(t) \rangle = 0 \qquad \langle z(t)z(s) \rangle = 0 \qquad \langle z(t)z^\ast(s) \rangle = \alpha(t-s)
+
+Here :math:`\langle \cdot \rangle` denotes the ensemble average.
+
+
+So far two different approaches to sample such processes have been implemented.
+
+Karhunen-Loève Expansion (KLE)
+------------------------------
+
+The Karhunen-Loève expansion makes use of the spectral representation of the auto correlation function viewed as a
+semi-positive integral kernel which allows for the decomposition.
+
+.. math:: \alpha(t-s) = \sum_n \lambda_n u_n(t) u_n^\ast(s)
+
+The Eigenvalues :math:`\lambda_n` and the Eigenfunction :math:`u_n(t)` of :math:`\alpha(t-s)` are defined in
+terms of the homogeneous Fredholm equation.
+
+.. math:: \int_0^T \mathrm{d} s \; \alpha(t-s) u_n(s) = \lambda_n u_n(t)
+
+For complex valued Gaussian distributed and independent random variables :math:`Y_n` with
+:math:`\langle Y_n \rangle = 0`,  :math:`\langle Y_n Y_m \rangle = 0` and
+:math:`\langle Y_n Y^\ast_m\rangle = \delta_{nm}` a stochastic process defined as
+
+.. math:: z(t) = \sum_n Y_n \sqrt{\lambda_n} u_n(t)
+
+obeys the required statistic. Note, the upper integral limit :math:`T` sets the time interval for which the
+stochastic process :math:`z(t) \; t \in [0,T]` is defined.
+
+The KLE approach is implemented by the class :py:class:`stocproc.StocProc_KLE`.
+It is numerically feasible if :math:`T` is not too large in comparision to a typical decay time of the
+auto correlation function.
+Implementation details can be found in the class documentation of :py:class:`stocproc.StocProc_KLE`.
+
+::
+
+    import stocproc as sp
+    import numpy as np
+    stp = sp.StocProc_KLE(alpha = lambda t: np.exp(- np.abs(t) - 1j*5*t),
+                          t_max = 2)
+    stp.new_process()
+    zt = stp()
+
+
+Fourier Transform Methods
+-------------------------
+
+This approach utilizes the relation of the auto correlation function and its Fourier transform.
+
+.. math:: \alpha(\tau) = \frac{1}{\pi} \int_{-\infty}^{\infty} \mathrm{d}\omega \; J(\omega) e^{-\mathrm{i}\omega\tau} \qquad J(\omega) \geq 0
+
+Discretizing the integral yields an approximate expression for the auto correlation function.
+
+.. math:: \alpha(\tau) \approx \sum_n w_n \frac{J(\omega_n)}{\pi} e^{-\mathrm{i}\omega_n\tau}
+
+For complex valued Gaussian distributed and independent random variables :math:`Y_n` with
+:math:`\langle Y_n \rangle = 0`,  :math:`\langle Y_n Y_m \rangle = 0` and
+:math:`\langle Y_n Y^\ast_m\rangle = \delta_{nm}` a stochastic process defined as
+
+.. math:: z(t) = \sum_n Y_n \sqrt{\frac{w_n J(\omega_n)}{\pi}} e^{-\mathrm{i}\omega_n t}
+
+obeys the required statistics up to an accuracy of the integral discretization.
+
+Fast Fourier Transform (FFT)
+````````````````````````````
+
+Equally distributed nodes :math:`\omega_n` allow for an evaluation of the stochastic process
+using the Fast Fourier Transform algorithm (see :py:class:`stocproc.StocProc_FFT` for an implementation).
+
+TanhSinh Intgeration (TanhSinh)
+```````````````````````````````
+
+For spectral densities :math:`J(\omega)` with a singularity at :math:`\omega=0` the TanhSinh integration
+scheme is more suitable. Such an implementation and its details can be found at :py:class:`stocproc.StocProc_TanhSinh`.
\ No newline at end of file
diff --git a/docs/source/stocproc.rst b/docs/source/stocproc.rst
deleted file mode 100644
index 33d9ffa..0000000
--- a/docs/source/stocproc.rst
+++ /dev/null
@@ -1 +0,0 @@
-.. automodule:: stocproc.stocproc
\ No newline at end of file
diff --git a/stocproc/__init__.py b/stocproc/__init__.py
index 7eacadf..a192b4a 100644
--- a/stocproc/__init__.py
+++ b/stocproc/__init__.py
@@ -1,150 +1,20 @@
-r"""
-Stochastic Process Module
-=========================
+__MAJOR__ = 1
+__MINOR__ = 0
+__PATCH__ = 0
 
+def version():
+    """semantic version string with format 'MAJOR.MINOR' (https://semver.org/)"""
+    return "{}.{}".format(__MAJOR__, __MINOR__)
 
-This module allows to sample Gaussian stochastic processes which are wide-sense stationary,
-continuous in time and complex valued. In particular for a given auto correlation function
-:math:`\alpha(\tau)` the stochastic process :math:`z(t)` obeys.
-
-.. math:: \langle z(t) \rangle = 0 \qquad \langle z(t)z(s) \rangle = 0 \qquad \langle z(t)z^\ast(s) \rangle = \alpha(t-s)
-
-Here :math:`\langle \cdot \rangle` denotes the ensemble average.
-
-
-So far two different approaches to sample such processes have been implemented.
-
-:doc:`Karhunen-Loève expansion (KLE) </StocProc_KLE>`
------------------------------------------------------
-
-The Karhunen-Loève expansion makes use of the spectral representation of the auto correlation function viewed as a
-semi-positive integral kernel which allows for the decomposition.
-
-.. math:: \alpha(t-s) = \sum_n \lambda_n u_n(t) u_n^\ast(s)
-
-The Eigenvalues :math:`\lambda_n` and the Eigenfunction :math:`u_n(t)` of :math:`\alpha(t-s)` are defined in
-terms of the homogeneous Fredholm equation.
-
-.. math:: \int_0^T \mathrm{d} s \; \alpha(t-s) u_n(s) = \lambda_n u_n(t)
-
-For complex valued Gaussian distributed and independent random variables :math:`Z_n` with
-:math:`\langle Z_n \rangle = 0`,  :math:`\langle Z_n Z_m \rangle = 0` and
-:math:`\langle Z_n Z^\ast_m\rangle = \delta_{nm}` a stochastic process defined as
-
-.. math:: z(t) = \sum_n Z_n \sqrt{\lambda_n} u_n(t)
-
-obeys the required statistic. Note, the upper integral limit :math:`T` sets the time interval for which the
-stochastic process :math:`z(t) \; t \in [0,T]` is defined.
-
-The KLE approach is implemented by the class :py:class:`stocproc.StocProc_KLE`.
-It is numerically feasible if :math:`T` is not too large in comparision to a typical decay time of the
-auto correlation function.
-Implementation details can be found in the class documentation of :py:class:`stocproc.StocProc_KLE`.
-
-::
-
-    import stocproc as sp
-    import numpy as np
-    stp = sp.StocProc_KLE(alpha = lambda t: np.exp(- np.abs(t) - 1j*5*t),
-                          t_max = 2)
-    stp.new_process()
-    zt = stp()
-
-
-:doc:`Fourier transform (FT) </StocProc_FFT>`
----------------------------------------------
-
-This approach utilizes the relation of the auto correlation function and its Fourier transform.
-
-.. math:: \alpha(\tau) = \frac{1}{\pi} \int_{-\infty}^{\infty} \mathrm{d}\omega \; J(\omega) e^{-\mathrm{i}\omega\tau} \qquad J(\omega) \geq 0
-
-Discretizing the integral yields an approximate expression for the auto correlation function.
-
-.. math:: \alpha(\tau) \approx \sum_n w_n \frac{J(\omega_n)}{\pi} e^{-\mathrm{i}\omega_n\tau}
-
-For complex valued Gaussian distributed and independent random variables :math:`Z_n` with
-:math:`\langle Z_n \rangle = 0`,  :math:`\langle Z_n Z_m \rangle = 0` and
-:math:`\langle Z_n Z^\ast_m\rangle = \delta_{nm}` a stochastic process defined as
-
-.. math:: z(t) = \sum_n Z_n \sqrt{\frac{w_n J(\omega_n)}{\pi}} e^{-\mathrm{i}\omega_n t}
-
-obeys the required statistics up to an accuracy of the integral discretization.
-
-Equally distributed nodes :math:`\omega_n` allow for an evaluation of the stochastic process
-using the Fast Fourier Transform algorithm (see :py:class:`stocproc.StocProc_FFT` for an implementation).
-
-For spectral densities :math:`J(\omega)` with a singularity at :math:`\omega=0` the TanhSinh integration
-scheme is more suitable. Such an implementation and its details can be found at :py:class:`stocproc.StocProc_TanhSinh`.
-
-
-
-To implement additional methods
-
-
-
-This module contains two different implementation for generating stochastic processes for a
-given auto correlation function (:doc:`Karhunen-Loève expansion </StocProc_KLE>`
-and :doc:`Fast-Fourier method </StocProc_FFT>`).
-Both methods are based on a time discrete process, however cubic
-spline interpolation is assured to be valid within a given tolerance.
-
-Documentation Overview
-----------------------
-
-.. toctree::
-   :maxdepth: 3
-
-   stocproc
-   example
-
-a Simple Example
-----------------
-
-The example will setup a process generator for an exponential auto correlation function
-and sample a single realization. ::
-
-    def lsd(w):
-        # Lorenzian spectral density
-        return 1/(1 + (w - _WC_)**2)
-
-
-    def exp_ac(t):
-        # exponential auto correlation function
-        # note there is a factor of one over pi in the
-        # definition of the auto correlation function:
-        # exp_ac(t) = 1/pi int_{-infty}^infty d w  lsd(w) exp(-i w t)
-        return np.exp(- np.abs(t) - 1j*_WC_*t)
-
-    _WC_ = 5
-    t_max = 10
-    print("setup process generator")
-    stp = sp.StocProc_FFT(spectral_density = lsd,
-                          t_max = t_max,
-                          bcf_ref = exp_ac,
-                          intgr_tol=1e-2,
-                          intpl_tol=1e-2):
-
-    print("generate single process")
-    stp.new_process()
-    zt = stp()    # get discrete process
-
-The full example can be found :doc:`here </example>`.
-
-.. image:: ../../examples/proc.*
-
-Averaging over 5000 samples yields the auto correlation function which is in good agreement
-with the exact auto correlation.
-
-.. image:: ../../examples/ac.*
-"""
-
-__version__ = "0.2.1"
+def version_full():
+    """semantic version string with format 'MAJOR.MINOR.PATCH' (https://semver.org/)"""
+    return "{}.{}.{}".format(__MAJOR__, __MINOR__, __PATCH__)
 
 import sys
 if sys.version_info.major < 3:
     raise SystemError("no support for Python 2")
 
-
+from .stocproc import loggin_setup
 from .stocproc import StocProc_FFT
 from .stocproc import StocProc_KLE
 from .stocproc import StocProc_TanhSinh
diff --git a/stocproc/method_fft.py b/stocproc/method_ft.py
similarity index 100%
rename from stocproc/method_fft.py
rename to stocproc/method_ft.py
diff --git a/stocproc/stocproc.py b/stocproc/stocproc.py
index 4100792..a9ad3cf 100644
--- a/stocproc/stocproc.py
+++ b/stocproc/stocproc.py
@@ -1,93 +1,86 @@
-"""
-Stochastic Process Generators
-=============================
-
-
-Karhunen-Loève expansion
-------------------------
-
-This method samples stochastic processes using Karhunen-Loève expansion and
-is implemented in the class :doc:`StocProc_KLE </StocProc_KLE>`.
-
-Setting up the class involves solving an eigenvalue problem which grows with
-the time interval the process is simulated on. Further generating a new process
-involves a multiplication with that matrix, therefore it scales quadratically with the
-time interval. Nonetheless it turns out that this method requires less random numbers
-than the Fast-Fourier method.
-
-
-Fast-Fourier method
--------------------
-
-In the class :doc:`StocProc_FFT </StocProc_FFT>` a method based on Fast-Fourier transform is
-used to sample stochastic processes.
-
-Setting up this class is quite efficient as it only calculates values of the
-associated spectral density. The number scales linear with the time interval of interest. However to achieve
-sufficient accuracy many of these values are required. As the generation of a new process is based on
-a Fast-Fouried-Transform over these values, this part is comparably lengthy.
-
-Abstract Base Class for a Stochastic Process Generator
-------------------------------------------------------
-
-:doc:`_absStocProc </StocProc_ABC>`.
-
-"""
 import abc
 import numpy as np
 import time
 
 from . import method_kle
-from . import method_fft
+from . import method_ft
 import fcSpline
 
 import logging
 log = logging.getLogger(__name__)
 
-ch = logging.StreamHandler()
-ch.setFormatter(logging.Formatter('%(name)s - %(levelname)s - %(message)s'))
+sh = logging.StreamHandler()
+sh.setFormatter(logging.Formatter('%(name)s - %(levelname)s - %(message)s'))
 
-def show_log(level = logging.INFO):
-    ch.setLevel(level)
+def loggin_setup(sh_level = logging.INFO, sp_log_level = logging.INFO,
+                 kle_log_level = logging.INFO, ft_log_level = logging.INFO):
+    """
+    controls the logging levels
 
-    log.addHandler(ch)
-    log.setLevel(level)
+    On loading the ``stocproc`` module this function is called with its default argument
+    which ammount to informative logging.
 
-    method_kle.log.addHandler(ch)
-    method_kle.log.setLevel(level)
+    :param sh_level: the logging level for the (output) StreamHandler
+    :param sp_log_level: the logging level of the stocproc module log
+    :param kle_log_level: the logging level of the KLE helper functions log (see :doc:`method_kle`)
+    :param ft_log_level: the logging level of the FT helper functions log (see :doc:`method_ft`)
+    """
+    sh.setLevel(sh_level)
 
-    method_fft.log.addHandler(ch)
-    method_fft.log.setLevel(level)
+    log.addHandler(sh)
+    log.setLevel(sp_log_level)
 
+    method_kle.log.addHandler(sh)
+    method_kle.log.setLevel(kle_log_level)
+
+    method_ft.log.addHandler(sh)
+    method_ft.log.setLevel(ft_log_level)
+
+loggin_setup()
 
 class _abcStocProc(abc.ABC):
     r"""
-    Abstract base class to stochastic process interface
-    
-    general work flow:
-        - Specify the time axis of interest [0, t_max] and it resolution (number of grid points), :math:`t_i = i \frac{t_max}{N_t-1}`.
-        - To evaluate the stochastic process at these points, a mapping from :math:`N_z` normal distributed 
-          random complex numbers with :math:`\langle y_i y_j^\ast \rangle = 2 \delta_{ij}`
-          to the stochastic process :math:`z_{t_i}` is needed and depends on the implemented method (:py:func:`calc_z`).
-        - A new process should be generated by calling :py:func:`new_process`.
-        - When the __call__ method is invoked the results will be interpolated between the :math:`z_{t_i}`.
-      
-    """
-    def __init__(self, t_max=None, num_grid_points=None, seed=None, t_axis=None, scale=1):
-        r"""
-            :param t_max: specify time axis as [0, t_max], if None, the times must be explicitly
-                given by t_axis
-            :param num_grid_points: number of equidistant times on that axis
-            :param seed: if not ``None`` set seed to ``seed``
-            :param t_axis: an explicit definition of times t_k (may be non equidistant)
-        """
-        if t_max is not None:
-            self.t = np.linspace(0, t_max, num_grid_points)
-        else:
-            self.t = t_axis
+    Interface definition for stochastic process implementations
 
-        self.t_max = self.t[-1]
-        self.num_grid_points = len(self.t)
+    A new implementation for a stochastic process generator should subclass :py:class:`_abcStocProc` and
+    overwrite :py:func:`calc_z` and :py:func:`get_num_y`.
+
+    Depending on the equally spaced times :math:`t_n = n \frac{t_{max}}{N-1}` with :math:`n = 0 \dots N-1`
+    (:math:`N` = number of grid points),
+    the function :py:func:`calc_z` should map :math:`M` independent complex valued and Gaussian
+    distributed random variables :math:`Y_m` (with :math:`\langle Y_m \rangle = 0 = \langle Y_m Y_{m'} \rangle` and
+    :math:`\langle Y_m Y^\ast_{m'} = \delta_{m m'}\rangle`) to the discrete time stochastic process :math:`z_n = z(t_n)`.
+
+    :py:func:`get_num_y` needs to return the number :math:`M` of random variables :math:`Y_m` required as
+    input for :py:func:`calc_z`.
+
+    Having implemented :py:func:`calc_z` and :py:func:`get_num_y` the :py:class:`_abcStocProc` provides
+    convenient functions such as:
+
+        - :py:func:`__call__`: evaluate the stochastic process for any time within the interval
+          :math:`[0, t_{max}]` using cubic spline interpolation
+        - :py:func:`get_time`: returns the times :math:`t_n`
+        - :py:func:`get_z`: returns the discrete stochastic process :math:`z_n`
+        - :py:func:`new_process`: draw new samples :math:`Y_m` and update :math:`z_n` as well as the
+          cubic spline interpolator
+        - :py:func:`set_scale`: set a scalar pre factor for the auto correlation function
+          which scales the stochastic process such that :math:`\langle z(t) z^\ast(s)\rangle = \text{scale} \; \alpha(t-s)`.
+
+    :param t_max: specifies the upper bound of the time interval
+    :param num_grid_points: number of grid points :math:`N`
+    :param seed: if not ``None`` seed the random number generator with ``seed``
+    :param t_axis: an explicit definition of times t_k (may be non equidistant)
+    :param scale: passes ``scale`` to :py:func:`set_scale`
+
+    Note: :py:func:`new_process` is **not** called on init. If you want to evaluate a particular
+    realization of the stocastic process, a new sample needs to be drawn by calling :py:func:`new_process`.
+    Otherwise a ``RuntimeError`` is raised.
+
+    """
+    def __init__(self, t_max=None, num_grid_points=None, seed=None, scale=1):
+        self.t_max = t_max
+        self.num_grid_points = num_grid_points
+        self.t = np.linspace(0, t_max, num_grid_points)
 
         self._z = None
         self._interpolator = None
@@ -101,13 +94,12 @@ class _abcStocProc(abc.ABC):
         log.debug("init StocProc with t_max {} and {} grid points".format(t_max, num_grid_points))
 
     def __call__(self, t=None):
-        r"""evaluates the stochastic process via spline interpolation of the discrete process :math:`z_k`
+        r"""Evaluates the stochastic process.
 
-        :param t: time to evaluate the stochastic process at, float of array of floats, if t is None
-            return the discrete process :math:`z_k` which corresponds to the times :math:`t_k` given by the
-            integration weights method
-        :return: a single complex value or a complex array of the shape of t that holds the values of
-            stochastic process at times t
+        If ``t`` is not ``None`` cubic spline interpolation is used to evaluate :math:`z(t)` based on the
+        discrete realization of the process :math:`z_n = z(t_n)`. ``t`` may be a single value or array like.
+
+        If ``t`` is ``None`` the discrete process :math:`z_n` is returned.
         """
         if self._z is None:
             raise RuntimeError("StocProc has NO random data, call 'new_process' to generate a new random process")
@@ -119,10 +111,13 @@ class _abcStocProc(abc.ABC):
 
     @abc.abstractmethod
     def calc_z(self, y):
-        r"""
-        maps the normal distributed complex valued random variables y to the stochastic process
-        
-        :return: the stochastic process, array of complex numbers 
+        r"""*abstract method*
+
+        An implementation needs to map :math:`M` independent complex valued and Gaussian
+        distributed random variables :math:`Y_m` (with :math:`\langle Y_m \rangle = 0 = \langle Y_m Y_{m'} \rangle` and
+        :math:`\langle Y_m Y^\ast_{m'} = \delta_{m m'}\rangle`) to the discrete time stochastic process :math:`z_n = z(t_n)`.
+
+        :return: the discrete time stochastic process :math:`z_n`, array of complex numbers
         """
         pass
     
@@ -132,41 +127,46 @@ class _abcStocProc(abc.ABC):
 
     @abc.abstractmethod
     def get_num_y(self):
-        r"""
-        :return: number of complex random variables needed to calculate the stochastic process 
+        r"""*abstract method*
+
+        An implementation needs to return the number :math:`M` of random variables :math:`Y_m` required as
+        input for :py:func:`calc_z`.
         """
         pass        
     
     def get_time(self):
-        r"""Returns the time :math:`t_k` corresponding to the values :math:`z_k`
-
-        These times are determined by the integration weights method.
-        """
+        r"""Returns the times :math:`t_n` for which the discrete time stochastic process may be evaluated."""
         return self.t
     
     def get_z(self):
-        r"""Returns the discrete process :math:`z_k`."""
+        r"""Returns the discrete time stochastic process :math:`z_n = z(t_n)`."""
         return self._z
     
     def new_process(self, y=None, seed=None):
-        r"""generate a new process by evaluating :py:func:`calc_z` with new random variables :math:`Y_i`
+        r"""Generate a new realization of the stochastic process.
 
-        :param y: independent normal distributed complex valued random variables with :math:`\sigma_{ij}^2 = \langle y_i y_j^\ast \rangle = 2 \delta_{ij}`
-        :param seed: if not None set seed to seed before generating samples
-        
-        When y is given use these random numbers as input for :py:func:`calc_z`
-        otherwise generate a new set of random numbers.
+        If ``seed`` is not ``None`` seed the random number generator with ``seed`` before drawing new random numbers.
+
+        If ``y`` is ``None`` draw new random numbers to generate the new realization. Otherwise use ``y`` as input to
+        generate the new realization.
         """
         t0 = time.time()
-        self._interpolator = None
+
+        # clean up old data
+        del self._interpolator
+        del self._z
+
         self._proc_cnt += 1
         if seed != None:
-            log.info("use fixed seed ({})for new process".format(seed))
+            log.info("use fixed seed ({}) for new process".format(seed))
             np.random.seed(seed)
         if y is None:
             #random complex normal samples
             y = np.random.normal(scale=self._one_over_sqrt_2, size = 2*self.get_num_y()).view(np.complex)
-        del self._z
+        else:
+            if len(y) != self.get_num_y():
+                raise RuntimeError("the length of 'y' ({}) needs to be {}".format(len(y), self.get_num_y()))
+
         self._z = self._calc_scaled_z(y)
         log.debug("proc_cnt:{} new process generated [{:.2e}s]".format(self._proc_cnt, time.time() - t0))
         t0 = time.time()
@@ -174,6 +174,10 @@ class _abcStocProc(abc.ABC):
         log.debug("created interpolator [{:.2e}s]".format(time.time() - t0))
         
     def set_scale(self, scale):
+        r"""
+        Set a scalar pre factor for the auto correlation function
+        which scales the stochastic process such that :math:`\langle z(t) z^\ast(s)\rangle = \text{scale} \; \alpha(t-s)`.
+        """
         self.scale = scale
         self.sqrt_scale = np.sqrt(scale)
 
@@ -294,18 +298,19 @@ class StocProc_KLE(_abcStocProc):
 
 
 class StocProc_FFT(_abcStocProc):
-    r"""Simulate Stochastic Process using FFT method
+    r"""Generate Stochastic Processes using the Fast Fourier Transform (FFT) method
 
-    This method uses the relation of the auto correlation to the non negative real valued
-    spectral density :math:`J(\omega)`. The integral can be approximated by a discrete integration scheme
+    This method uses the relation of the auto correlation function ``alpha`` to the non negative real valued
+    spectral density (``spectral_density``) :math:`J(\omega)`.
+    The integral can be approximated by a discrete integration scheme
 
     .. math::
-        \alpha(\tau) = \int_{\omega_\mathrm{min}}^{\omega_\mathrm{max}} \mathrm{d}\omega \, \frac{J(\omega)}{\pi} e^{-\mathrm{i}\omega \tau}
+        \alpha(\tau) = \int_{-\infty}^{\infty} \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} \omega_k \tau}
 
     where the weights :math:`\omega_k` depend on the particular integration scheme. For a process defined as
 
-    .. math:: Z(t) = \sum_{k=0}^{N-1} \sqrt{\frac{w_k J(\omega_k)}{\pi}} Y_k \exp^{-\mathrm{i}\omega_k t}
+    .. math:: 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 :math:`Y_k` such that :math:`\langle Y_k \rangle = 0`,
     :math:`\langle Y_k Y_{k'}\rangle = 0` and :math:`\langle Y_k Y^\ast_{k'}\rangle = \delta_{k,k'}`
@@ -313,90 +318,97 @@ class StocProc_FFT(_abcStocProc):
 
     .. math::
         \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)}
+            \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}
 
-    To calculate :math:`Z(t)` the Discrete Fourier Transform (DFT) can be applied as follows:
+    To calculate :math:`z(t)` the Discrete Fourier Transform (DFT) can be applied as follows:
 
-    .. math:: 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}
+    .. math:: 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}
+
+    However, this requires that :math:`\omega_k` takes the form :math:`\omega_k = \omega_\mathrm{min} + k \Delta \omega`
+    with :math:`\Delta \omega = (\omega_\mathrm{max} - \omega_\mathrm{min}) / (N-1)` which limits
+    the integration schemes to those with equidistant nodes.
 
-    Here :math:`\omega_k` has to take the form :math:`\omega_k = \omega_\mathrm{min} + k \Delta \omega` and
-    :math:`\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 :math:`\Delta t` has to be chosen such that
     :math:`2\pi = N \Delta \omega \Delta t` holds.
-    Since :math:`J(\omega)` is real it follows that :math:`X(t_l) = X^\ast(t_{N-l})`.
+    Since :math:`J(\omega)` is real it follows that :math:`z(t_l) = z^\ast(t_{N-l})`.
     For that reason the stochastic process has only :math:`(N+1)/2` (odd :math:`N`) and
     :math:`(N/2 + 1)` (even :math:`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].
+    To generate a process with given auto correlation function on the interval :math:`[0, t_{max}]`
+    requires that the auto correlation function approximation is valid for all :math:`t` in :math:`[0, t_{max}]`.
 
     This is ensured by automatically determining the number of sumands N and the integral
     boundaries :math:`\omega_\mathrm{min}` and :math:`\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 :math:`t_l \in [0, t_\mathrm{max}]`.
+    within the tolerance intgr_tol for all discrete :math:`t_l \in [0, t_{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 :math:`\Delta t`
+    due to the interpolation is controlled by the parameter ``intpl_tol``. The maximum time step :math:`\Delta t`
     is chosen such that the interpolated valued at each half step :math:`t_i + \Delta t /2` differs at
-    most intpl_tol from the exact value of the auto correlation function.
+    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 :math:`\Delta t`
     criterion from the interpolation is met.
 
+    See :py:func:`stocproc.method_ft.calc_ab_N_dx_dt` for implementation details on how the
+    tolerance criterion are met. Since the pre calculation may become time consuming the :py:class:`StocProc_FFT`
+    class can be pickled and unpickled. To identify a particular instance a unique key is formed by the tuple
+    ``(alpha, t_max, intgr_tol, intpl_tol)``.
+    It is advisable to use :py:func:`get_key` with keyword arguments to generate such a tuple.
+
 
     :param spectral_density: the spectral density :math:`J(\omega)` as callable function object
     :param t_max: :math:`[0,t_\mathrm{max}]` is the interval for which the process will be calculated
-    :param bcf_ref: a callable which evaluates the Fourier integral exactly
+    :param alpha: a callable which evaluates the Fourier integral exactly
     :param intgr_tol: tolerance for the integral approximation
     :param intpl_tol: tolerance for the interpolation
     :param seed: if not None, use this seed to seed the random number generator
     :param negative_frequencies: if False, keep :math:`\omega_\mathrm{min} = 0` otherwise
        find a negative :math:`\omega_\mathrm{min}` appropriately just like :math:`\omega_\mathrm{max}`
 
-    .. todo::
-       implement bcf_ref = None and use numeric integration as default
-
-
     """
-    def __init__(self, spectral_density, t_max, bcf_ref, intgr_tol=1e-2, intpl_tol=1e-2,
+    def __init__(self, spectral_density, t_max, alpha, intgr_tol=1e-2, intpl_tol=1e-2,
                  seed=None, negative_frequencies=False, scale=1):
-        self.key = bcf_ref, t_max, intgr_tol, intpl_tol
+        self.key = self.get_key(alpha=alpha, t_max=t_max, intgr_tol=intgr_tol, intpl_tol=intpl_tol)
                 
         if not negative_frequencies: 
             log.info("non neg freq only")
-            a, b, N, dx, dt = method_fft.calc_ab_N_dx_dt(integrand = spectral_density,
-                                                         intgr_tol = intgr_tol,
-                                                         intpl_tol = intpl_tol,
-                                                         t_max     = t_max,
-                                                         ft_ref    = lambda tau:bcf_ref(tau)*np.pi,
-                                                         opt_b_only= True)
+            a, b, N, dx, dt = method_ft.calc_ab_N_dx_dt(integrand = spectral_density,
+                                                        intgr_tol = intgr_tol,
+                                                        intpl_tol = intpl_tol,
+                                                        t_max     = t_max,
+                                                        ft_ref    = lambda tau:alpha(tau)*np.pi,
+                                                        opt_b_only= True)
         else:
             log.info("use neg freq")
-            a, b, N, dx, dt = method_fft.calc_ab_N_dx_dt(integrand = spectral_density,
-                                                         intgr_tol = intgr_tol,
-                                                         intpl_tol = intpl_tol,
-                                                         t_max     = t_max,
-                                                         ft_ref    = lambda tau:bcf_ref(tau)*np.pi,
-                                                         opt_b_only= False)
+            a, b, N, dx, dt = method_ft.calc_ab_N_dx_dt(integrand = spectral_density,
+                                                        intgr_tol = intgr_tol,
+                                                        intpl_tol = intpl_tol,
+                                                        t_max     = t_max,
+                                                        ft_ref    = lambda tau:alpha(tau)*np.pi,
+                                                        opt_b_only= False)
 
-        assert abs(2*np.pi - N*dx*dt) < 1e-12
+        d = abs(2*np.pi - N*dx*dt)
+        if d >= 1e-12:
+            log.fatal('method_ft.calc_ab_N_dx_dt returned inconsistent data!')
+            raise RuntimeError('d = {:.3e} < 1e-12 FAILED!'.format(d))
 
-        print("Fourier Integral Boundaries: [{:.3e}, {:.3e}]".format(a,b))
-        print("Number of Nodes            : {}".format(N))
-        print("yields dx                  : {:.3e}".format(dx))
-        print("yields dt                  : {:.3e}".format(dt))
-        print("yields t_max               : {:.3e}".format( (N-1)*dt))
+        log.info("Fourier Integral Boundaries: [{:.3e}, {:.3e}]".format(a,b))
+        log.info("Number of Nodes            : {}".format(N))
+        log.info("yields dx                  : {:.3e}".format(dx))
+        log.info("yields dt                  : {:.3e}".format(dt))
+        log.info("yields t_max               : {:.3e}".format( (N-1)*dt))
 
         num_grid_points = int(np.ceil(t_max/dt))+1
 
-        assert num_grid_points <= N
+        if num_grid_points > N:
+            log.fatal("num_grid_points and number of points used for FFT inconsistent!")
+            raise RuntimeError("num_grid_points = {} <= N_DFT = {}  FAILED!".format(num_grid_points, N))
 
         t_max = (num_grid_points-1)*dt
-        
         super().__init__(t_max           = t_max, 
                          num_grid_points = num_grid_points, 
                          seed            = seed,
@@ -407,8 +419,12 @@ class StocProc_FFT(_abcStocProc):
         self.omega_min_correction = np.exp(-1j*(a+dx/2)*self.t)   #self.t is from the parent class
 
     @staticmethod
-    def get_key(t_max, bcf_ref, intgr_tol=1e-2, intpl_tol=1e-2):
-        return bcf_ref, t_max, intgr_tol, intpl_tol
+    def get_key(t_max, alpha, intgr_tol=1e-2, intpl_tol=1e-2):
+        """
+        Returns the tuple ``(alpha, t_max, intgr_tol, intpl_tol)`` which uniquely identifies a particular
+        :py:class:`StocProc_FFT` instance
+        """
+        return alpha, t_max, intgr_tol, intpl_tol
 
     def __getstate__(self):
         return self.yl, self.num_grid_points, self.omega_min_correction, self.t_max, self._seed, self.scale, self.key
@@ -421,20 +437,20 @@ class StocProc_FFT(_abcStocProc):
                          scale           = scale)
             
     def calc_z(self, y):
-        r"""calculate
+        r"""Calculate the discrete time stochastic process using FFT algorithm
 
         .. math::
-            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)
+            z_n = z(t_n) = e^{-\mathrm{i}\omega_\mathrm{min} t_n} \mathrm{FFT}\left( \sqrt{\frac{w_k J(\omega_k)}{\pi}} Y_k \right)
 
-        and return values with :math:`t_l < t_\mathrm{max}`
+        and return values :math:`z_n` with :math:`t_n <= t_\mathrm{max}`.
         """
         z_fft = np.fft.fft(self.yl * y)
         z = z_fft[0:self.num_grid_points] * self.omega_min_correction
         return z
 
     def get_num_y(self):
-        r"""The number of independent random variables Y is given by the number of discrete times
-        :math:`t_l < t_\mathrm{max}` from the Fourier Transform
+        r"""The number of independent random variables :math:`Y_m` is given by the number of discrete nodes
+        used by the Fast Fourier Transform algorithm. 
         """
         return len(self.yl)
 
@@ -452,13 +468,13 @@ class StocProc_TanhSinh(_abcStocProc):
             log.info("get_dt_for_accurate_interpolation, please wait ...")
             try:
                 ft_ref = lambda tau: bcf_ref(tau) * np.pi
-                c = method_fft.find_integral_boundary(lambda tau: np.abs(ft_ref(tau)) / np.abs(ft_ref(0)),
+                c = method_ft.find_integral_boundary(lambda tau: np.abs(ft_ref(tau)) / np.abs(ft_ref(0)),
                                                       intgr_tol, 1, 1e6, 0.777)
             except RuntimeError:
                 c = t_max
 
             c = min(c, t_max)
-            dt_tol = method_fft.get_dt_for_accurate_interpolation(t_max=c,
+            dt_tol = method_ft.get_dt_for_accurate_interpolation(t_max=c,
                                                                   tol=intpl_tol,
                                                                   ft_ref=ft_ref)
             log.info("requires dt < {:.3e}".format(dt_tol))
@@ -468,7 +484,7 @@ class StocProc_TanhSinh(_abcStocProc):
         N = int(np.ceil(t_max/dt_tol))+1
         log.info("yields N = {}".format(N))
 
-        wmax = method_fft.find_integral_boundary(spectral_density, tol=intgr_tol/4, ref_val=1, max_val=1e6, x0=0.777)
+        wmax = method_ft.find_integral_boundary(spectral_density, tol=intgr_tol/4, ref_val=1, max_val=1e6, x0=0.777)
         log.info("wmax:{}".format(wmax))
 
         tau = np.linspace(0, t_max, 35)
@@ -476,7 +492,7 @@ class StocProc_TanhSinh(_abcStocProc):
         d = intgr_tol+1
         while d > intgr_tol:
             n *= 2
-            I = method_fft.fourier_integral_TanhSinh(f = lambda w: spectral_density(w)/np.pi,
+            I = method_ft.fourier_integral_TanhSinh(f = lambda w: spectral_density(w)/np.pi,
                                                      x_max=wmax,
                                                      n=n,
                                                      tau_l=tau)
@@ -487,7 +503,7 @@ class StocProc_TanhSinh(_abcStocProc):
 
         tau = np.linspace(0, (N-1)*dt_tol, N)
         log.info("perform numeric check of entire time axis [{},{}] n:{}".format(0, (N-1)*dt_tol, N))
-        num_FT = method_fft.fourier_integral_TanhSinh(f = lambda w: spectral_density(w)/np.pi,
+        num_FT = method_ft.fourier_integral_TanhSinh(f = lambda w: spectral_density(w)/np.pi,
                                                       x_max=wmax,
                                                       n=n,
                                                       tau_l=tau)
@@ -511,7 +527,7 @@ class StocProc_TanhSinh(_abcStocProc):
             plt.show()
 
         assert d <= intgr_tol, "d:{}, intgr_tol:{}".format(d, intgr_tol)
-        yk, wk = method_fft.get_x_w_and_dt(n, wmax)
+        yk, wk = method_ft.get_x_w_and_dt(n, wmax)
         self.omega_k = yk
         self.fl = np.sqrt(wk*spectral_density(self.omega_k)/np.pi)