.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_data_analysis/estimate_stochastic_processes/plot_estimate_spectral_density_function.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        Click :ref:`here <sphx_glr_download_auto_data_analysis_estimate_stochastic_processes_plot_estimate_spectral_density_function.py>`
        to download the full example code

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_data_analysis_estimate_stochastic_processes_plot_estimate_spectral_density_function.py:


Estimate a spectral density function
====================================

.. GENERATED FROM PYTHON SOURCE LINES 6-19

The objective of this example is to estimate the spectral density
function :math:`S` from data, which can be a sample of time series or one time
series.

The following example illustrates the case where the available data is a
sample of :math:`10^3` realizations of the process, defined on the time grid
:math:`[0, 102.3]`, discretized every :math:`\Delta t = 0.1`. The spectral model of
the process is the Cauchy model parameterized by :math:`\underline{\lambda}=(5)`
and :math:`\underline{a}=(3)`.

The figure draws the graph of the real spectral
model and its estimation from the sample of time series.


.. GENERATED FROM PYTHON SOURCE LINES 21-26

.. code-block:: default

    import openturns as ot
    import openturns.viewer as viewer
    from matplotlib import pylab as plt
    ot.Log.Show(ot.Log.NONE)








.. GENERATED FROM PYTHON SOURCE LINES 27-28

generate some data

.. GENERATED FROM PYTHON SOURCE LINES 28-47

.. code-block:: default


    # Create the time grid
    # In the context of the spectral estimate or Fourier transform use,
    # we use data blocs with size of form 2^p
    tMin = 0.
    tstep = 0.1
    size = 2**12
    tgrid = ot.RegularGrid(tMin, tstep, size)

    # We fix the parameter of the Cauchy model
    amplitude = [5.0]
    scale = [3.0]
    model = ot.CauchyModel(amplitude, scale)
    process = ot.SpectralGaussianProcess(model, tgrid)

    # Get a time series or a sample of time series
    tseries = process.getRealization()
    sample = process.getSample(1000)








.. GENERATED FROM PYTHON SOURCE LINES 48-49

Build a spectral model factory

.. GENERATED FROM PYTHON SOURCE LINES 49-53

.. code-block:: default

    segmentNumber = 10
    overlapSize = 0.3
    factory = ot.WelchFactory(ot.Hann(), segmentNumber, overlapSize)








.. GENERATED FROM PYTHON SOURCE LINES 54-55

Estimation on a TimeSeries or on a ProcessSample

.. GENERATED FROM PYTHON SOURCE LINES 55-58

.. code-block:: default

    estimatedModel_TS = factory.build(tseries)
    estimatedModel_PS = factory.build(sample)








.. GENERATED FROM PYTHON SOURCE LINES 59-60

Change the filtering window

.. GENERATED FROM PYTHON SOURCE LINES 60-62

.. code-block:: default

    factory.setFilteringWindows(ot.Hamming())








.. GENERATED FROM PYTHON SOURCE LINES 63-64

Get the frequencyGrid

.. GENERATED FROM PYTHON SOURCE LINES 64-67

.. code-block:: default

    frequencyGrid = ot.SpectralGaussianProcess(
        estimatedModel_PS, tgrid).getFrequencyGrid()








.. GENERATED FROM PYTHON SOURCE LINES 68-98

.. code-block:: default


    # With the model, we want to compare values
    # We compare values computed with theoritical values
    plotSample = ot.Sample(frequencyGrid.getN(), 3)

    # Loop of comparison ==> data are saved in plotSample
    for k in range(frequencyGrid.getN()):
        freq = frequencyGrid.getStart() + k * frequencyGrid.getStep()
        plotSample[k, 0] = freq
        plotSample[k, 1] = abs(estimatedModel_PS(freq)[0, 0])
        plotSample[k, 2] = abs(model(freq)[0, 0])

    # Some cosmetics : labels, legend position, ...
    graph = ot.Graph("Estimated spectral function - Validation", "Frequency",
                     "Spectral density function", True, "topright", 1.0, ot.GraphImplementation.LOGY)

    # The first curve is the estimate density as function of frequency
    curve1 = ot.Curve(plotSample.getMarginal([0, 1]))
    curve1.setColor('blue')
    curve1.setLegend('estimate model')

    # The second curve is the theoritical density as function of frequency
    curve2 = ot.Curve(plotSample.getMarginal([0, 2]))
    curve2.setColor('red')
    curve2.setLegend('Cauchy model')

    graph.add(curve1)
    graph.add(curve2)
    view = viewer.View(graph)
    plt.show()



.. image-sg:: /auto_data_analysis/estimate_stochastic_processes/images/sphx_glr_plot_estimate_spectral_density_function_001.png
   :alt: Estimated spectral function - Validation
   :srcset: /auto_data_analysis/estimate_stochastic_processes/images/sphx_glr_plot_estimate_spectral_density_function_001.png
   :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** ( 0 minutes  1.486 seconds)


.. _sphx_glr_download_auto_data_analysis_estimate_stochastic_processes_plot_estimate_spectral_density_function.py:


.. only :: html

 .. container:: sphx-glr-footer
    :class: sphx-glr-footer-example



  .. container:: sphx-glr-download sphx-glr-download-python

     :download:`Download Python source code: plot_estimate_spectral_density_function.py <plot_estimate_spectral_density_function.py>`



  .. container:: sphx-glr-download sphx-glr-download-jupyter

     :download:`Download Jupyter notebook: plot_estimate_spectral_density_function.ipynb <plot_estimate_spectral_density_function.ipynb>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_