Example 2: ECG signal denoising with the SASS algorithm

This example shows the use of the sparsity-assisted signal smoothing (SASS) algorithm for ECG filtering. In this example, the QRS waveform is modeled as piecewise quadratic, so we use K = 3 in SASS.

Ivan Selesnick
selesi@poly.edu
2013

Contents

Start

clc
clear

printme = @(filename) print('-dpdf', sprintf('figures/Example2_%s', filename));

Make ECG signal

Use 'ecgsyn.m' by P. E. McSharry, G. D. Clifford, L. Tarassenko, and L. A. Smith. A dynamical model for generating synthetic electrocardiogram signals. Trans. on Biomed. Eng., 50(3):289-294, March 2003. http://www.physionet.org/physiotools/ecgsyn/

addpath ecg
randn('state', 0);
rand('state', 1);               % Set rand/randn states for reproducibility

Fs = 256;                       % Fs : sampling frequency (samples/sec)

ecg_10 = ecgsyn(Fs, 10);        % Simulate 10 heart beats
ECG sampled at 256 Hz
Approximate number of heart beats: 10
Measurement noise amplitude: 0 
Heart rate mean: 60 bpm
Heart rate std: 1 bpm
LF/HF ratio: 0.5
Internal sampling frequency: 512
      P  Q  R  S  T
ti = [-1.22173 -0.261799 0 0.261799 1.74533] radians
ai = [1.2 -5 30 -7.5 0.75]
bi = [0.25 0.1 0.1 0.1 0.4]
Integrating dynamical system

d =

     4

M = 2;                          % M : number of beats
ecg = ecg_10(Fs/2+(1:M*Fs));    % Extract M beats from simulated ECG

N = length(ecg);
n = 0:N-1;

sigma = 0.1;                    % sigma : noise standard deviation

noise = sigma * randn(N, 1);    % noise : white Gaussian noise

data = ecg + noise;             % data : noisy ECG

Display ECG signals

ax1 = [0 N -1.0 2.0];

figure(1)
clf
subplot(2,1,1)
plot(n, ecg)
title('ECG (noise-free)')
axis(ax1)

subplot(2,1,2)
plot(n, data)
title('ECG + noise');
axis(ax1)

Perform preprocessing

The SASS algorithm may introduce transients at the beginning and end of the signal. To reduce the transients, replace first and last samples of the data by a low-order polynomial approximation.

r = 2;
M = 15;
y = preproc(r, M, data);

figure(1)
clf
subplot(2, 1, 1)
plot(n, y)
axis(ax1)
title('Preprocessed data')

Filter matrices

Define filter matrices for SASS.

d = 2;                          % filter is of order 2d
fc = 0.03;                      % fc : cut-off frequency (0 < fc < 0.5) (cycles/sample)
K = 3;                          % K : order of sparse derivative

% Banded filter matrices:
[A, B, B1, D, a, b, b1, H1norm HTH1norm] = ABfilt(d, fc, N, K);

H = @(x) [nan(d,1); A\(B*x); nan(d,1)];     % H: high-pass filter
L = @(x) x - H(x);                          % L: low-pass filter

txt_filt = sprintf('d = %d, fc = %.3f', d, fc);

Low-pass filter

x_lpf = L(y);

txt_LPF = sprintf('Low-pass filtering (%s)', txt_filt);

err = ecg - x_lpf;
RMSE_LPF = sqrt(mean(err(K+1:end-K).^2))

figure(1)
clf
subplot(2, 1, 1)
plot(n, x_lpf)
title(txt_LPF)
axis(ax1)
RMSE_LPF =

    0.1172

SASS (L1 norm)

beta = 3;
lam = beta * sigma * HTH1norm;  % lam : regularization parameter

Nit = 100;                      % Nit : number of iterations

[x_L1, cost_L1, u_L1, v_L1] = sass(y, d, fc, K, lam, 'L1', [], Nit);
% x : denoised signal
% u : sparse order-K derivative
% v : inv(A)*B1*u
% cost : cost function history

txt_L1 = sprintf('SASS (L1, K = %d, \\lambda = %.2f, %s)', K, lam, txt_filt)

err = ecg - x_L1;
RMSE_L1 = sqrt(mean(err(K+1:end-K).^2))
0 incorrectly locked zeros detected.

txt_L1 =

SASS (L1, K = 3, \lambda = 2.74, d = 2, fc = 0.030)


RMSE_L1 =

    0.0382

Cost function

figure(1)
clf
plot(cost_L1)
title('Cost function history');
xlabel('Iteration')

Denoised signal

figure(1)
clf
subplot(2, 1, 1)
plot(n, x_L1)
title(txt_L1)
axis(ax1)
axnote(sprintf('RMSE = %.4f', RMSE_L1))

Optimality scatter plot

figure(1)
clf
plot_scatter(A, B, B1, y, u_L1, lam, 'L1', []);
title(sprintf('L1 optimality scatter plot (%d iterations)', Nit))
printme('L1_scatter')

Components

figure(1)
clf
plot_components(n, data, x_lpf, x_L1, u_L1, v_L1, txt_L1);
printme('L1_components')

SASS (atan penalty)

The arctangent (atan) penalty function is non-convex

% rho : Controls non-convexity of penalty (0 <= rho <= 1)
rho = 0.5;

% initialize with L1 solution
[x_atan1, cost_atan1, u_atan1, v_atan1, a] = sass(y, d, fc, K, lam, 'atan', rho, Nit, u_L1);
3 incorrectly locked zeros detected.

Correct for zero-locking

[x_atan, cost_atan, u_atan, v_atan, a] = sass2(y, d, fc, K, lam, 'atan', rho, Nit, u_L1);

cost_delta = cost_atan1(end) - cost_atan(end);

fprintf('Zero-locking correction improved cost function value by %g\n', cost_delta);
Run 1: 3 incorrectly locked zeros detected.
Run 2: 0 incorrectly locked zeros detected.
Zero-locking correction improved cost function value by 0.000277121

Display scatter plot before and after zero-locking correction

figure(2)
clf
subplot(2, 1, 1)
plot_scatter(A, B, B1, y, u_atan1, lam, 'atan', a);
ax_scatter = [-0.005 0.005 0.5 1.3];
axis(ax_scatter)
title('Without zero-locking correction (not optimal)')

subplot(2, 1, 2)
plot_scatter(A, B, B1, y, u_atan, lam, 'atan', a);
axis(ax_scatter)
title('With zero-locking correction (locally optimal)')

orient tall
printme('atan_scatter_detail')
txt_atan = sprintf('SASS (atan, K = %d, \\lambda = %.2f, \\rho = %.2f, %s)', K, lam, rho, txt_filt);

err = ecg - x_atan;
RMSE_atan = sqrt(mean(err(K+1:end-K).^2))
RMSE_atan =

    0.0306

Cost function

figure(1)
clf
plot(cost_atan)
title('Cost function history');
xlabel('Iteration')

Denoised signal

figure(1)
clf
subplot(2, 1, 1)
plot(n, x_atan)
title(txt_atan)
axis(ax1)
axnote(sprintf('RMSE = %.4f', RMSE_atan))

Optimality scatter plot

figure(1)
clf
plot_scatter(A, B, B1, y, u_atan, lam, 'atan', a);
title(sprintf('atan optimality scatter plot (%d iterations)', Nit))
printme('atan_scatter')

Components

figure(1)
clf
plot_components(n, data, x_lpf, x_atan, u_atan, v_atan, txt_atan);
printme('atan_components')

Save figure to pdf file

ax1 = [0 N -0.5 1.5];

set(0, 'DefaultAxesFontSize', 8);

figure(2)
clf

subplot(5, 1, 1)
plot(n, ecg)
title('ECG (noise-free)');
axis(ax1)

subplot(5, 1, 2)
plot(n, data)
title('ECG + noise');
axis(ax1)

subplot(5, 1, 3)
plot(n, x_lpf);
title(txt_LPF);
axis(ax1)
axnote(sprintf('RMSE = %.4f', RMSE_LPF))

subplot(5, 1, 4)
plot(n, x_L1);
title(txt_L1);
axis(ax1)
axnote(sprintf('RMSE = %.4f', RMSE_L1))

subplot(5, 1, 5)
plot(n, x_atan);
title(txt_atan);
axis(ax1)
axnote(sprintf('RMSE = %.4f', RMSE_atan))

orient tall
print -dpdf figures/Example2_fig1


set(0, 'DefaultAxesFontSize', 'remove');