{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n\n# ECG Data Compressive Sensing\n :depth: 2\n :local:\n\nIn this example, we demonstrate the compressive sensing of ECG data\nand reconstruction using Block Sparse Bayesian Learning (BSBL).\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Configure JAX to work with 64-bit floating point precision. \nfrom jax.config import config\nconfig.update(\"jax_enable_x64\", True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's import necessary libraries\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import timeit\nimport jax \nimport numpy as np\nimport jax.numpy as jnp\n# CR-Suite libraries\nimport cr.nimble as crn\nimport cr.nimble.dsp as crdsp\nimport cr.sparse.dict as crdict\nimport cr.sparse.plots as crplot\nimport cr.sparse.block.bsbl as bsbl\n\n# Sample data\nfrom scipy.misc import electrocardiogram\n# Plotting\nimport matplotlib.pyplot as plt\n# Miscellaneous\nfrom scipy.signal import detrend, butter, filtfilt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Test signal\nSciPy includes a test electrocardiogram signal\nwhich is a 5 minute long electrocardiogram (ECG), \na medical recording of the electrical activity of the heart, \nsampled at 360 Hz.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ecg = electrocardiogram()\n# Sampling frequency in Hz\nfs = 360\n# We shall only process a part of the signal in this demo\nN = 400\nx = ecg[:N]\nt = np.arange(N) * (1/fs)\nfig, ax = plt.subplots(figsize=(16,4))\nax.plot(t, x);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Preprocessing\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Remove the linear trend from the signal\nx = detrend(x)\n## bandpass filter\n# lower cutoff frequency\nf1 = 5\n# upper cutoff frequency\nf2 = 40\n# passband in normalized frequency\nWn = np.array([f1, f2]) * 2 / fs\n# butterworth filter\nfn = 3\nfb, fa = butter(fn, Wn, 'bandpass')\nx = filtfilt(fb,fa,x)\nfig, ax = plt.subplots(figsize=(16,4))\nax.plot(t, x);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compressive Sensing at 70%\nWe choose the compression ratio (M/N) to be 0.7\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "CR = 0.70\nM = int(N * CR)\nprint(f'M={M}, N={N}, CR={CR}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Sensing matrix\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "Phi = crdict.gaussian_mtx(crn.KEY0, M, N)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Measurements\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "y = Phi @ x\nfig, ax = plt.subplots(figsize=(16, 4))\nax.plot(y);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Sparse Recovery with BSBL\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "options = bsbl.bsbl_bo_options(y, max_iters=20)\nstart = timeit.default_timer()\nsol = bsbl.bsbl_bo_np_jit(Phi, y, 25, options=options)\nstop = timeit.default_timer()\nprint(f'Reconstruction time: {stop - start:.2f} sec', )\nprint(sol)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Recovered signal\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "x_hat = sol.x\nprint(f'SNR: {crn.signal_noise_ratio(x, x_hat):.2f} dB, PRD: {crn.percent_rms_diff(x, x_hat):.1f}%')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the original and recovered signals\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ax = crplot.h_plots(2)\nax[0].plot(x)\nax[1].plot(x_hat)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compressive Sensing at 50%\nLet us now increase the compression\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "CR = 0.50\nM = int(N * CR)\nprint(f'M={M}, N={N}, CR={CR}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Sensing matrix\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "Phi = crdict.gaussian_mtx(crn.KEY0, M, N)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Measurements\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "y = Phi @ x\nfig, ax = plt.subplots(figsize=(16, 4))\nax.plot(y);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Sparse Recovery with BSBL\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "options = bsbl.bsbl_bo_options(y, max_iters=20)\nstart = timeit.default_timer()\nsol = bsbl.bsbl_bo_np_jit(Phi, y, 25, options=options)\nstop = timeit.default_timer()\nprint(f'Reconstruction time: {stop - start:.2f} sec', )\nprint(sol)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Recovered signal\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "x_hat = sol.x\nprint(f'SNR: {crn.signal_noise_ratio(x, x_hat):.2f} dB, PRD: {crn.percent_rms_diff(x, x_hat):.1f}%')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the original and recovered signals\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ax = crplot.h_plots(2)\nax[0].plot(x)\nax[1].plot(x_hat)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.13" } }, "nbformat": 4, "nbformat_minor": 0 }