.. 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

        :ref:`Go to the end <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-27

.. code-block:: Python

    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 28-29

generate some data

.. GENERATED FROM PYTHON SOURCE LINES 29-48

.. code-block:: Python


    # 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.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 49-50

Build a spectral model factory

.. GENERATED FROM PYTHON SOURCE LINES 50-54

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 55-56

Estimation on a TimeSeries or on a ProcessSample

.. GENERATED FROM PYTHON SOURCE LINES 56-59

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 60-61

Change the filtering window

.. GENERATED FROM PYTHON SOURCE LINES 61-63

.. code-block:: Python

    factory.setFilteringWindows(ot.Hamming())








.. GENERATED FROM PYTHON SOURCE LINES 64-65

Get the frequencyGrid

.. GENERATED FROM PYTHON SOURCE LINES 65-67

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 68-105

.. code-block:: Python


    # 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,
        "upper right",
        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






.. _sphx_glr_download_auto_data_analysis_estimate_stochastic_processes_plot_estimate_spectral_density_function.py:

.. only:: html

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

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

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

    .. 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>`