.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/ecg_compression.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_ecg_compression.py: ECG Data Compression ===================================== .. contents:: :depth: 2 :local: In this example, we demonstrate the compression of ECG data using biorthogonal wavelets. We shall first go through different steps of the wavelet compression algorithm one by one. We will then show a unified encoder decoder function. This example is adapted from the sample code provided in `nerajbobra/wavelet-based-ecg-compression `_. Do refer to its documentation to get a general sense of the compression algorithm. This implementation is significantly different and optimized. .. GENERATED FROM PYTHON SOURCE LINES 23-29 .. code-block:: default # Configure JAX to work with 64-bit floating point precision. from jax.config import config config.update("jax_enable_x64", True) .. GENERATED FROM PYTHON SOURCE LINES 30-31 Let's import necessary libraries .. GENERATED FROM PYTHON SOURCE LINES 31-48 .. code-block:: default import jax import numpy as np import jax.numpy as jnp # CR-Suite libraries import cr.nimble as crn import cr.nimble.dsp as crdsp import cr.wavelets as crwt from cr.wavelets.codec import * from cr.wavelets.plots import * # Sample data from scipy.misc import electrocardiogram # Plotting import matplotlib.pyplot as plt # Miscellaneous from scipy.signal import detrend .. GENERATED FROM PYTHON SOURCE LINES 49-51 Encoder Configuration ------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 51-63 .. code-block:: default # Number of samples per signal NUM_SAMPLES_BLOCK = 3000 # Name of wavelet to be used for signal decomposition WAVELET_NAME = 'bior4.4' # Number of levels of decomposition WAVELET_LEVEL = 5 # Fraction of energy to be preserved in each decomposition level WAVELET_ENERGY_THRESHOLDS = [0.999, 0.97, 0.85, 0.85, 0.85, 0.85] # Maximum percentage root mean square difference acceptable in quantization MAX_PRD = 40 # percent .. GENERATED FROM PYTHON SOURCE LINES 64-70 Test signal ------------------------------ SciPy includes a test electrocardiogram signal which is a 5 minute long electrocardiogram (ECG), a medical recording of the electrical activity of the heart, sampled at 360 Hz. .. GENERATED FROM PYTHON SOURCE LINES 70-80 .. code-block:: default ecg = electrocardiogram() # Sampling frequency in Hz fs = 360 # We shall only process one signal in this demo signal = ecg[:NUM_SAMPLES_BLOCK] n = len(signal) t = np.arange(n) * (1/fs) fig, ax = plt.subplots(figsize=(16,4)) ax.plot(t, signal); .. image-sg:: /gallery/images/sphx_glr_ecg_compression_001.png :alt: ecg compression :srcset: /gallery/images/sphx_glr_ecg_compression_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/cr-wavelets/checkouts/stable/examples/ecg_compression.py:70: DeprecationWarning: scipy.misc.electrocardiogram has been deprecated in SciPy v1.10.0; and will be completely removed in SciPy v1.12.0. Dataset methods have moved into the scipy.datasets module. Use scipy.datasets.electrocardiogram instead. ecg = electrocardiogram() [] .. GENERATED FROM PYTHON SOURCE LINES 81-82 Preprocessing .. GENERATED FROM PYTHON SOURCE LINES 82-88 .. code-block:: default # Remove the linear trend from the signal signal = detrend(signal) fig, ax = plt.subplots(figsize=(16,4)) ax.plot(t, signal); .. image-sg:: /gallery/images/sphx_glr_ecg_compression_002.png :alt: ecg compression :srcset: /gallery/images/sphx_glr_ecg_compression_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 89-91 Encoder ------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 93-94 Let us perform a wavelet decomposition of this signal .. GENERATED FROM PYTHON SOURCE LINES 94-100 .. code-block:: default wavelet = crwt.to_wavelet(WAVELET_NAME) coeffs = crwt.wavedec(signal, wavelet, level=WAVELET_LEVEL) fig, axes = plt.subplots(6, 1, figsize=(16,25)) plot_decomposition(fig, axes, coeffs) .. image-sg:: /gallery/images/sphx_glr_ecg_compression_003.png :alt: cA5, cD5, cD4, cD3, cD2, cD1 :srcset: /gallery/images/sphx_glr_ecg_compression_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 101-102 Let us check the quality of the reconstruction .. GENERATED FROM PYTHON SOURCE LINES 102-106 .. code-block:: default rec = crwt.waverec(coeffs, wavelet) fig, ax = plt.subplots(figsize=(16,5)) plot_2_signals(fig, ax, signal, rec, fs=fs) .. image-sg:: /gallery/images/sphx_glr_ecg_compression_004.png :alt: ecg compression :srcset: /gallery/images/sphx_glr_ecg_compression_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 107-108 Percentage root mean square difference .. GENERATED FROM PYTHON SOURCE LINES 108-110 .. code-block:: default print(crn.prd(signal, rec)) .. rst-class:: sphx-glr-script-out .. code-block:: none 1.0158245596205641e-05 .. GENERATED FROM PYTHON SOURCE LINES 111-114 Let us threshold different levels of decomposition by keeping as few large coefficients as required to meet the target energy fraction (separately for each level) .. GENERATED FROM PYTHON SOURCE LINES 114-119 .. code-block:: default th_coeffs, binmaps = threshold(coeffs, WAVELET_ENERGY_THRESHOLDS) # Let us plot them to see the coefficients which have been zeroed out fig, axes = plt.subplots(6, 1, figsize=(16,25)) plot_2_decompositions(fig, axes, coeffs, th_coeffs, label2='Thresholded') .. image-sg:: /gallery/images/sphx_glr_ecg_compression_005.png :alt: cA5, cD5, cD4, cD3, cD2, cD1 :srcset: /gallery/images/sphx_glr_ecg_compression_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 120-121 Let us remove the zero entries from thresholded coefficients .. GENERATED FROM PYTHON SOURCE LINES 121-126 .. code-block:: default nz_th_coeffs = remove_zeros(th_coeffs, binmaps) # number of non-zero coefficients left at each level for c in nz_th_coeffs: print(len(c)) .. rst-class:: sphx-glr-script-out .. code-block:: none 95 46 15 24 149 484 .. GENERATED FROM PYTHON SOURCE LINES 127-128 Let us check how much loss we have suffered due to thresholding .. GENERATED FROM PYTHON SOURCE LINES 128-139 .. code-block:: default # Add back the zeros tmp = add_zeros(nz_th_coeffs, binmaps) # Perform reconstruction using thresholded coefficients th_rec = crwt.waverec(tmp, wavelet) fig, ax = plt.subplots(figsize=(16,5)) plot_2_signals(fig, ax, signal, th_rec, fs=fs) # Compute the percentage root mean square difference print(crn.prd(signal, th_rec)) .. image-sg:: /gallery/images/sphx_glr_ecg_compression_006.png :alt: ecg compression :srcset: /gallery/images/sphx_glr_ecg_compression_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none 16.568432073643812 .. GENERATED FROM PYTHON SOURCE LINES 140-141 Let us now scale the nonzero thresholded coefficients to 0-1 range .. GENERATED FROM PYTHON SOURCE LINES 141-146 .. code-block:: default scaled_coeffs, shifts, scales = scale_to_0_1(nz_th_coeffs) # Let us plot to verify that they are indeed now in 0-1 range fig, axes = plt.subplots(6, 1, figsize=(16,25)) plot_decomposition(fig, axes, scaled_coeffs, stem=True) .. image-sg:: /gallery/images/sphx_glr_ecg_compression_007.png :alt: cA5, cD5, cD4, cD3, cD2, cD1 :srcset: /gallery/images/sphx_glr_ecg_compression_007.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 147-149 Let us quantize the scaled coefficients so that we meet the maximum PRD criteria. .. GENERATED FROM PYTHON SOURCE LINES 149-159 .. code-block:: default (quantized_coeffs, num_bits, cur_prd) = quantize_to_prd_target(signal, WAVELET_NAME, WAVELET_LEVEL, scaled_coeffs, shifts, scales, binmaps, MAX_PRD) # Let us look at the quantized coefficients at each decomposition level fig, axes = plt.subplots(6, 1, figsize=(16,25)) plot_decomposition(fig, axes, quantized_coeffs, stem=True) # Check how many bits are being used per sample and the achieved PRD print(num_bits, cur_prd) .. image-sg:: /gallery/images/sphx_glr_ecg_compression_008.png :alt: cA5, cD5, cD4, cD3, cD2, cD1 :srcset: /gallery/images/sphx_glr_ecg_compression_008.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none 4 18.17883530994624 .. GENERATED FROM PYTHON SOURCE LINES 160-161 Merge the coefficients and binary maps of different levels to single arrays .. GENERATED FROM PYTHON SOURCE LINES 161-166 .. code-block:: default combined_coeffs = combine_arrays(quantized_coeffs) combined_binmaps = combine_arrays(binmaps) print(len(combined_coeffs), len(combined_binmaps)) .. rst-class:: sphx-glr-script-out .. code-block:: none 813 3041 .. GENERATED FROM PYTHON SOURCE LINES 167-177 We are now ready with all the data that needs to be transmitted to the decoder. This includes: - The quantized nonzero coefficients - The binary maps - The shifts and scales for each level used during scaling - The number of bits per sample used for quantization We shall use the coding function which will compress all of this data into a packed ``bitarray``. .. GENERATED FROM PYTHON SOURCE LINES 177-179 .. code-block:: default result = encode_cbss_to_bits(combined_coeffs, combined_binmaps, shifts, scales, num_bits) .. GENERATED FROM PYTHON SOURCE LINES 180-181 Check the number of bits used in compressing the whole block .. GENERATED FROM PYTHON SOURCE LINES 181-182 .. code-block:: default print(len(result)) .. rst-class:: sphx-glr-script-out .. code-block:: none 7108 .. GENERATED FROM PYTHON SOURCE LINES 183-184 We can also check the average number of bits per sample .. GENERATED FROM PYTHON SOURCE LINES 184-185 .. code-block:: default print(len(result)/NUM_SAMPLES_BLOCK) .. rst-class:: sphx-glr-script-out .. code-block:: none 2.3693333333333335 .. GENERATED FROM PYTHON SOURCE LINES 186-188 The MIT-BIH database is encoded at 11-bits per sample. We can use this to compute the compression ratio .. GENERATED FROM PYTHON SOURCE LINES 188-191 .. code-block:: default compression_ratio = len(signal) * 11 / len(result) print(compression_ratio) .. rst-class:: sphx-glr-script-out .. code-block:: none 4.642656162070906 .. GENERATED FROM PYTHON SOURCE LINES 192-198 Decoding ------------------------------ We note that the decoder has access to only the encoded bitstream and the encoder configuration parameters. It doesn't have access to any other intermediate data that was produced during encoding. .. GENERATED FROM PYTHON SOURCE LINES 200-203 The first step is to extract all the data from the packed bitarray. Note that the number of bits per sample used for quantization is also being read from the bitstream. .. GENERATED FROM PYTHON SOURCE LINES 203-214 .. code-block:: default (dec_c_coeffs, dec_c_binmaps, dec_shifts, dec_scales, dec_qbits) = decode_cbss_from_bits( WAVELET_NAME, WAVELET_LEVEL, NUM_SAMPLES_BLOCK, result) # Since we have encoding steps data available, we can # cross check to see if the data extraction from the # bitstream happened correctly. print(np.allclose(shifts, dec_shifts)) print(np.allclose(scales, dec_scales)) print(np.allclose(combined_binmaps, dec_c_binmaps)) print(np.allclose(combined_coeffs, dec_c_coeffs)) .. rst-class:: sphx-glr-script-out .. code-block:: none True True True True .. GENERATED FROM PYTHON SOURCE LINES 215-217 Let us now split the coefficients and binary maps to different decomposition levels .. GENERATED FROM PYTHON SOURCE LINES 217-221 .. code-block:: default dec_coeffs, dec_binmaps = split_coefs_binmaps( WAVELET_NAME, WAVELET_LEVEL, NUM_SAMPLES_BLOCK, dec_c_coeffs, dec_c_binmaps) .. GENERATED FROM PYTHON SOURCE LINES 222-223 Perform inverse quantization of nonzero coefficients .. GENERATED FROM PYTHON SOURCE LINES 223-225 .. code-block:: default inv_quant_coeffs = inv_quantize_1(dec_coeffs, dec_qbits) .. GENERATED FROM PYTHON SOURCE LINES 226-227 Perform descaling of the coefficients from [0, 1] range to their original ranges .. GENERATED FROM PYTHON SOURCE LINES 227-229 .. code-block:: default dec_unscaled_coeffs = descale_from_0_1(inv_quant_coeffs, dec_shifts, dec_scales) .. GENERATED FROM PYTHON SOURCE LINES 230-231 Add back the zero entries using the binary maps .. GENERATED FROM PYTHON SOURCE LINES 231-232 .. code-block:: default dec_coeffs = add_zeros(dec_unscaled_coeffs, dec_binmaps) .. GENERATED FROM PYTHON SOURCE LINES 233-234 Perform reconstruction of the signal from the decoded wavelet decomposition .. GENERATED FROM PYTHON SOURCE LINES 234-238 .. code-block:: default dec_reconstructed = crwt.waverec(dec_coeffs, wavelet) # Plot the original signal and reconstructed signal together fig, ax = plt.subplots(figsize=(16,5)) plot_2_signals(fig, ax, signal, dec_reconstructed, fs=fs) .. image-sg:: /gallery/images/sphx_glr_ecg_compression_009.png :alt: ecg compression :srcset: /gallery/images/sphx_glr_ecg_compression_009.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 239-240 Measure the percentage root mean square difference .. GENERATED FROM PYTHON SOURCE LINES 240-244 .. code-block:: default print(crn.prd(signal, dec_reconstructed)) .. rst-class:: sphx-glr-script-out .. code-block:: none 18.178835029200155 .. GENERATED FROM PYTHON SOURCE LINES 245-252 Full CODEC ------------------------------ Carefully carrying out individual steps in the compression and reconstruction is hard. It would be great if the library can provide simple functions which wrap all the encoding and decoding operations together. .. GENERATED FROM PYTHON SOURCE LINES 254-257 We can use the ``build_codec`` function to build an encoder and decoder function based on the encoder configuration parameters .. GENERATED FROM PYTHON SOURCE LINES 257-260 .. code-block:: default encoder, decoder = build_codec( WAVELET_NAME, WAVELET_LEVEL, NUM_SAMPLES_BLOCK, MAX_PRD, WAVELET_ENERGY_THRESHOLDS) .. GENERATED FROM PYTHON SOURCE LINES 261-262 Let us use the encoder to compress the signal .. GENERATED FROM PYTHON SOURCE LINES 262-263 .. code-block:: default bits = encoder(signal) .. GENERATED FROM PYTHON SOURCE LINES 264-266 Check that the encoder function did exactly the same thing as our step by step procedure above. .. GENERATED FROM PYTHON SOURCE LINES 266-267 .. code-block:: default print(result == bits) .. rst-class:: sphx-glr-script-out .. code-block:: none True .. GENERATED FROM PYTHON SOURCE LINES 268-269 Now reconstruct the signal from the encoded bitstream .. GENERATED FROM PYTHON SOURCE LINES 269-270 .. code-block:: default signal_rec = decoder(bits) .. GENERATED FROM PYTHON SOURCE LINES 271-272 Plot the original signal against the decoded signal .. GENERATED FROM PYTHON SOURCE LINES 272-274 .. code-block:: default fig, ax = plt.subplots(figsize=(16,5)) plot_2_signals(fig, ax, signal, signal_rec, fs=fs) .. image-sg:: /gallery/images/sphx_glr_ecg_compression_010.png :alt: ecg compression :srcset: /gallery/images/sphx_glr_ecg_compression_010.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 275-276 Measure the percentage root mean square difference .. GENERATED FROM PYTHON SOURCE LINES 276-277 .. code-block:: default print(crn.prd(signal, signal_rec)) .. rst-class:: sphx-glr-script-out .. code-block:: none 18.178835029200155 .. GENERATED FROM PYTHON SOURCE LINES 278-279 Compute the compression ratio .. GENERATED FROM PYTHON SOURCE LINES 279-284 .. code-block:: default compression_ratio = len(signal) * 11 / len(bits) print(compression_ratio) .. rst-class:: sphx-glr-script-out .. code-block:: none 4.642656162070906 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 11.085 seconds) .. _sphx_glr_download_gallery_ecg_compression.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: ecg_compression.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: ecg_compression.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_