.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/06_denoise/06_01_1d_fractional_brownian_motion.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. or to run this example in your browser via JupyterLite or Binder .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_06_denoise_06_01_1d_fractional_brownian_motion.py: 1D Fractional Brownian Motion ============================= We use the smooth module to smooth a 1D fractional Brownian motion signal. .. GENERATED FROM PYTHON SOURCE LINES 13-15 Imports ------- .. GENERATED FROM PYTHON SOURCE LINES 15-22 .. code-block:: Python import math import numpy as np import matplotlib.pyplot as plt from splineops.denoise.fBmper import fBmper from splineops.denoise.smoothing_spline import smoothing_spline .. GENERATED FROM PYTHON SOURCE LINES 23-28 1D Fractional Brownian Motion ----------------------------- A realization of an fBm process of length N is generated and corrupted with noise. The sequence is then denoised and oversampled by a factor m using the optimal fractional-spline estimator. .. GENERATED FROM PYTHON SOURCE LINES 28-94 .. code-block:: Python # Define program constants m = 4 # Upsampling factor N = 256 # Number of samples # Default values default_H = 0.7 default_SNRmeas = 20.0 default_verify = '0' # Enter Hurst parameter [0 < H < 1] (default: 0.7) > H = 0.7 # Enter measurement SNR at the mid-point (t = {N/2}) [dB] (default: 20.0) > SNRmeas = 20.0 # Create pseudo-fBm signal epsH = 1 t0, y0 = fBmper(epsH, H, m, N) Ch = epsH ** 2 / (math.gamma(2 * H + 1) * np.sin(np.pi * H)) POWmid = Ch * (N / 2) ** (2 * H) # Theoretical fBm variance at the midpoint # Measurement: downsample and add noise t = t0[::m] y = y0[::m] sigma = np.sqrt(POWmid) / (10 ** (SNRmeas / 20)) noise = np.random.randn(N) noise = sigma * noise / np.sqrt(np.mean(noise ** 2)) y_noisy = y + noise # Find smoothing spline fit lambda_ = (sigma / epsH) ** 2 gamma_ = H + 0.5 ts, ys = smoothing_spline(y_noisy, lambda_, m, gamma_) # Add non-stationary correction term cnn = np.concatenate(([1], np.zeros(N - 1))) # Normalized white noise autocorrelation tes, r = smoothing_spline(cnn, lambda_, m, gamma_) r = r * ys[0] / r[0] y_est = ys - r # Calculate MSE and SNR MSE0 = np.mean(noise ** 2) # Measurement MSE MSE = np.mean((y_est[::m] - y0[::m]) ** 2) # Denoised sequence MSE MSEm = np.mean((y_est - y0) ** 2) # Denoised and oversampled signal MSE SNR0 = 10 * np.log10(POWmid / MSE0) SNR = 10 * np.log10(POWmid / MSE) SNRm = 10 * np.log10(POWmid / MSEm) print(f'Number of measurements is {N}, oversampling factor is {m}.') print(f'mSNR (SNR at the mid-point) of the measured sequence is {SNR0:.2f} dB.') print(f'mSNR improvement of the denoised sequence is {SNR - SNR0:.2f} dB.') print(f'mSNR improvement of the denoised and oversampled sequence is {SNRm - SNR0:.2f} dB.') # Plot the results plt.figure() plt.plot(t0[:len(t0)//2], y0[:len(y0)//2], 'k', label='fBm') plt.plot(t[:len(t)//2], y_noisy[:len(y_noisy)//2], 'k+:', label='Noisy fBm samples') plt.plot(ts[:len(ts)//2], ys[:len(ys)//2], 'r--', label='Stationary estimation') plt.plot(tes[:len(tes)//2], y_est[:len(y_est)//2], 'r', label='Non-stationary estimation') plt.legend() plt.title(f'Estimation of fBm (H = {H}, ε_H^2 = {epsH}, σ_N^2 = {sigma ** 2:.4f})') plt.xlabel('Time') plt.ylabel('B_H') plt.tight_layout() plt.show() .. image-sg:: /auto_examples/06_denoise/images/sphx_glr_06_01_1d_fractional_brownian_motion_001.png :alt: Estimation of fBm (H = 0.7, ε_H^2 = 1, σ_N^2 = 8.8707) :srcset: /auto_examples/06_denoise/images/sphx_glr_06_01_1d_fractional_brownian_motion_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Number of measurements is 256, oversampling factor is 4. mSNR (SNR at the mid-point) of the measured sequence is 20.00 dB. mSNR improvement of the denoised sequence is 7.57 dB. mSNR improvement of the denoised and oversampled sequence is 7.71 dB. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.157 seconds) .. _sphx_glr_download_auto_examples_06_denoise_06_01_1d_fractional_brownian_motion.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/splineops/splineops.github.io/main?urlpath=lab/tree/notebooks_binder/auto_examples/06_denoise/06_01_1d_fractional_brownian_motion.ipynb :alt: Launch binder :width: 150 px .. container:: lite-badge .. image:: images/jupyterlite_badge_logo.svg :target: ../../lite/lab/index.html?path=auto_examples/06_denoise/06_01_1d_fractional_brownian_motion.ipynb :alt: Launch JupyterLite :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 06_01_1d_fractional_brownian_motion.ipynb <06_01_1d_fractional_brownian_motion.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 06_01_1d_fractional_brownian_motion.py <06_01_1d_fractional_brownian_motion.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 06_01_1d_fractional_brownian_motion.zip <06_01_1d_fractional_brownian_motion.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_