.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/0004.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 .. rst-class:: sphx-glr-example-title .. _sphx_glr_gallery_0004.py: .. _gallery:0004: Complex Sinusoids+Complex-Spikes in Dirac-Fourier Basis ========================================================== .. contents:: :depth: 2 :local: The Complex Sinusoids+Spikes signal in this example consists of a mixture of 3 different complex sinusoids with different amplitudes/phases and 120 different complex spikes (with normally distributed amplitudes for both real and imaginary parts and randomly chosen locations). The sinusoid part of the signal has a sparse representation in the Fourier basis. The Spikes part of the signal has a sparse representation in the Dirac(Identity/Standard) basis. Thus, the mixture Sinusoid+Spikes signal has a sparse representation (of 123 nonzero coefficients) in the Dirac-Fourier two ortho basis. Note that the spikes are normally distributed. Some of the spikes have extremely low amplitudes. They may be missed by a recovery algorithm depending on convergence thresholds. See also: * :ref:`api:problems` * :ref:`api:lop` .. GENERATED FROM PYTHON SOURCE LINES 37-48 .. code-block:: default # Configure JAX to work with 64-bit floating point precision. from jax.config import config config.update("jax_enable_x64", True) import jax.numpy as jnp import cr.nimble as crn import cr.sparse.plots as crplot .. GENERATED FROM PYTHON SOURCE LINES 49-53 Setup ------------------------------ We shall construct our test signal and dictionary using our test problems module. .. GENERATED FROM PYTHON SOURCE LINES 53-59 .. code-block:: default from cr.sparse import problems prob = problems.generate('complex:sinusoid-spikes:dirac-fourier', c=3, k=120) fig, ax = problems.plot(prob) .. image-sg:: /gallery/images/sphx_glr_0004_001.png :alt: Real part of the signal, Imaginary part of the signal, Real part of the coefficients, Imaginary part of the coefficients :srcset: /gallery/images/sphx_glr_0004_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 60-61 Let us access the relevant parts of our test problem .. GENERATED FROM PYTHON SOURCE LINES 61-70 .. code-block:: default # The sparsifying basis linear operator A = prob.A # The Complex Sinusoids+Spikes signal b0 = prob.b # The sparse representation of the signal in the dictionary x0 = prob.x .. GENERATED FROM PYTHON SOURCE LINES 71-73 Check how many coefficients in the sparse representation are sufficient to capture 99.9% of the energy of the signal .. GENERATED FROM PYTHON SOURCE LINES 73-75 .. code-block:: default print(crn.num_largest_coeffs_for_energy_percent(x0, 99.9)) .. rst-class:: sphx-glr-script-out .. code-block:: none 102 .. GENERATED FROM PYTHON SOURCE LINES 76-78 This number gives us an idea about the required sparsity to be configured for greedy pursuit algorithms. .. GENERATED FROM PYTHON SOURCE LINES 80-83 Sparse Recovery using Subspace Pursuit ------------------------------------------- We shall use subspace pursuit to reconstruct the signal. .. GENERATED FROM PYTHON SOURCE LINES 83-87 .. code-block:: default import cr.sparse.pursuit.sp as sp # We will first try to estimate a 100-sparse representation sol = sp.solve(A, b0, 100) .. GENERATED FROM PYTHON SOURCE LINES 88-89 This utility function helps us quickly analyze the quality of reconstruction .. GENERATED FROM PYTHON SOURCE LINES 89-90 .. code-block:: default problems.analyze_solution(prob, sol) .. rst-class:: sphx-glr-script-out .. code-block:: none m: 1024, n: 2048 b_norm: original: 71.089 reconstruction: 71.047 SNR: 29.32 dB x_norm: original: 71.188 reconstruction: 71.241 SNR: 29.29 dB Sparsity: original: 102, reconstructed: 93, overlap: 93, ratio: 0.912 Iterations: 20 .. GENERATED FROM PYTHON SOURCE LINES 91-93 It takes 20 iterations to converge and 76 of the largest 78 entries have been correctly identified. .. GENERATED FROM PYTHON SOURCE LINES 95-96 We will now try to estimate a 150-sparse representation .. GENERATED FROM PYTHON SOURCE LINES 96-98 .. code-block:: default sol = sp.solve(A, b0, 150) problems.analyze_solution(prob, sol) .. rst-class:: sphx-glr-script-out .. code-block:: none m: 1024, n: 2048 b_norm: original: 71.089 reconstruction: 71.089 SNR: 293.43 dB x_norm: original: 71.188 reconstruction: 71.188 SNR: 293.70 dB Sparsity: original: 102, reconstructed: 102, overlap: 102, ratio: 1.000 Iterations: 3 .. GENERATED FROM PYTHON SOURCE LINES 99-102 We have correctly detected all the 78 most significant entries Let us check if we correctly decoded all the nonzero entries in the sparse representation x .. GENERATED FROM PYTHON SOURCE LINES 102-104 .. code-block:: default problems.analyze_solution(prob, sol, perc=100) .. rst-class:: sphx-glr-script-out .. code-block:: none m: 1024, n: 2048 b_norm: original: 71.089 reconstruction: 71.089 SNR: 293.43 dB x_norm: original: 71.188 reconstruction: 71.188 SNR: 293.70 dB Sparsity: original: 123, reconstructed: 123, overlap: 123, ratio: 1.000 Iterations: 3 .. GENERATED FROM PYTHON SOURCE LINES 105-106 The estimated sparse representation .. GENERATED FROM PYTHON SOURCE LINES 106-107 .. code-block:: default x = sol.x .. GENERATED FROM PYTHON SOURCE LINES 108-109 Let us reconstruct the signal from this sparse representation .. GENERATED FROM PYTHON SOURCE LINES 109-111 .. code-block:: default b = prob.reconstruct(x) .. GENERATED FROM PYTHON SOURCE LINES 112-113 Let us visualize the original and reconstructed representation .. GENERATED FROM PYTHON SOURCE LINES 113-126 .. code-block:: default def plot_representations(x0, x): ax = crplot.h_plots(4, height=2) ax[0].stem(jnp.real(x0), markerfmt='.') ax[0].set_title('Original real part') ax[1].stem(jnp.real(x), markerfmt='.') ax[1].set_title('Reconstructed real part') ax[2].stem(jnp.imag(x0), markerfmt='.') ax[2].set_title('Original imaginary part') ax[3].stem(jnp.imag(x), markerfmt='.') ax[3].set_title('Reconstructed imaginary part') plot_representations(x0, x) .. image-sg:: /gallery/images/sphx_glr_0004_002.png :alt: Original real part, Reconstructed real part, Original imaginary part, Reconstructed imaginary part :srcset: /gallery/images/sphx_glr_0004_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 127-128 Let us visualize the original and reconstructed signal .. GENERATED FROM PYTHON SOURCE LINES 128-141 .. code-block:: default def plot_signals(b0, b): ax = crplot.h_plots(4, height=2) ax[0].plot(jnp.real(b0)) ax[0].set_title('Original real part') ax[1].plot(jnp.real(b)) ax[1].set_title('Reconstructed real part') ax[2].plot(jnp.imag(b0)) ax[2].set_title('Original imaginary part') ax[3].plot(jnp.imag(b)) ax[3].set_title('Reconstructed imaginary part') plot_signals(b0, b) .. image-sg:: /gallery/images/sphx_glr_0004_003.png :alt: Original real part, Reconstructed real part, Original imaginary part, Reconstructed imaginary part :srcset: /gallery/images/sphx_glr_0004_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 142-145 Sparse Recovery using Compressive Sampling Matching Pursuit --------------------------------------------------------------- We shall now use compressive sampling matching pursuit to reconstruct the signal. .. GENERATED FROM PYTHON SOURCE LINES 145-150 .. code-block:: default import cr.sparse.pursuit.cosamp as cosamp # We will try to estimate a 150-sparse representation sol = cosamp.solve(A, b0, 150) problems.analyze_solution(prob, sol, perc=100) .. rst-class:: sphx-glr-script-out .. code-block:: none m: 1024, n: 2048 b_norm: original: 71.089 reconstruction: 71.089 SNR: 292.93 dB x_norm: original: 71.188 reconstruction: 71.188 SNR: 292.70 dB Sparsity: original: 123, reconstructed: 123, overlap: 123, ratio: 1.000 Iterations: 4 .. GENERATED FROM PYTHON SOURCE LINES 151-152 The estimated sparse representation .. GENERATED FROM PYTHON SOURCE LINES 152-153 .. code-block:: default x = sol.x .. GENERATED FROM PYTHON SOURCE LINES 154-155 Let us reconstruct the signal from this sparse representation .. GENERATED FROM PYTHON SOURCE LINES 155-158 .. code-block:: default b = prob.reconstruct(x) .. GENERATED FROM PYTHON SOURCE LINES 159-160 Let us visualize the original and reconstructed representation .. GENERATED FROM PYTHON SOURCE LINES 160-163 .. code-block:: default plot_representations(x0, x) .. image-sg:: /gallery/images/sphx_glr_0004_004.png :alt: Original real part, Reconstructed real part, Original imaginary part, Reconstructed imaginary part :srcset: /gallery/images/sphx_glr_0004_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 164-165 Let us visualize the original and reconstructed signal .. GENERATED FROM PYTHON SOURCE LINES 165-168 .. code-block:: default plot_signals(b0, b) .. image-sg:: /gallery/images/sphx_glr_0004_005.png :alt: Original real part, Reconstructed real part, Original imaginary part, Reconstructed imaginary part :srcset: /gallery/images/sphx_glr_0004_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 169-171 Sparse Recovery using SPGL1 --------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 171-176 .. code-block:: default import cr.sparse.cvx.spgl1 as crspgl1 options = crspgl1.SPGL1Options() sol = crspgl1.solve_bp_jit(A, b0, options=options) problems.analyze_solution(prob, sol, perc=100) .. rst-class:: sphx-glr-script-out .. code-block:: none m: 1024, n: 2048 b_norm: original: 71.089 reconstruction: 71.089 SNR: 118.43 dB x_norm: original: 71.188 reconstruction: 71.188 SNR: 113.81 dB Sparsity: original: 123, reconstructed: 123, overlap: 123, ratio: 1.000 Iterations: 41 n_times: 44, n_trans: 42 .. GENERATED FROM PYTHON SOURCE LINES 177-178 The estimated sparse representation .. GENERATED FROM PYTHON SOURCE LINES 178-179 .. code-block:: default x = sol.x .. GENERATED FROM PYTHON SOURCE LINES 180-181 Let us reconstruct the signal from this sparse representation .. GENERATED FROM PYTHON SOURCE LINES 181-184 .. code-block:: default b = prob.reconstruct(x) .. GENERATED FROM PYTHON SOURCE LINES 185-186 Let us visualize the original and reconstructed representation .. GENERATED FROM PYTHON SOURCE LINES 186-189 .. code-block:: default plot_representations(x0, x) .. image-sg:: /gallery/images/sphx_glr_0004_006.png :alt: Original real part, Reconstructed real part, Original imaginary part, Reconstructed imaginary part :srcset: /gallery/images/sphx_glr_0004_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 190-191 Let us visualize the original and reconstructed signal .. GENERATED FROM PYTHON SOURCE LINES 191-194 .. code-block:: default plot_signals(b0, b) .. image-sg:: /gallery/images/sphx_glr_0004_007.png :alt: Original real part, Reconstructed real part, Original imaginary part, Reconstructed imaginary part :srcset: /gallery/images/sphx_glr_0004_007.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 195-204 Comments --------- * With K=100, SP recovery is slightly inaccurate. It also takes more iterations to converge. * Both SP and CoSaMP correctly recover the signal in 3 iterations if the sparsity is specified properly (150 > 123). * SPGL1 converges in 41 iterations but correctly discovers the support without requiring any input about expected sparsity. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 16.582 seconds) .. _sphx_glr_download_gallery_0004.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 0004.py <0004.py>` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 0004.ipynb <0004.ipynb>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_