Example 2: Local polynomial approximation + total variation filtering

Simultaneous local least-square polynomial approximation and total variation filtering

Ivan Selesnick,
Polytechnic Institute of NYU
December 2011
Email: selesi@poly.edu
Reference: Polynomial Smoothing of Time Series with Additive Step Discontinuities
I. W. Selesnick, S. Arnold, and V. R. Dantham

Contents

Start

clear
% close all

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

Create signals

Create simulated signal (smooth signal with additive step discontinuities)

N = 300;
n = (1:N)';

s1 = 2*(n < 0.3*N) + 1*(n > 0.6*N);         % step function
s2 = sin(0.021*pi*n);                       % smooth function
s = s1 + s2;                                % total signal

randn('state',0);                           % Initialize randn so that example can be exactly reproduced
sigma = 0.3;
noise = sigma*randn(N,1);
y = s + noise;                              % noisy signal

figure(2)
clf
subplot(2,1,1)
plot(y,'.k')
title('Noisy data');

printme('data')

LoPATV filtering

Apply local polynomial approximation + TV filtering (LoPATV filter) to noisy data

deg = 2;            % deg : degree of polynomial
P = 40;             % P : block overlap
L = 50;             % L : block length
lambda = 8;         % lambda : regularization parameter

(N-L)/(L-P)+1       % This is the number of blocks - it should be an
                    % integer, otherwise the data will be truncated

Nit = 500;          % Nit : number of iterations
mu0 = 10;           % mu0 : augmented Lagrangian parameter
mu = 1;             % mu : augmented Lagrangian parameter

[x, p, cost] = lopatv(y, L, P, deg, lambda, Nit, mu0, mu);

rmse_lopatv = sqrt(mean((s - x - p).^2));

fprintf('LoPATV: RMSE = %f\n', rmse_lopatv)
ans =

    26

LoPATV: RMSE = 0.086840

Display cost function history. The cost function flattens out as the algorithm converges.

figure(1)
clf
plot(cost)
title('Cost function history (LoPATV filtering algorithm)')
xlabel('Iteration')
printme('LoPATV_convergence')

Display TV component and TV-compensated data

figure(1)
clf
subplot(2,1,1)
plot(x, 'black')
title('Calculated TV component (LoPATV)');

subplot(2,1,2)
plot(y - x, 'k.')
title('TV-compensated data (LoPATV)');

printme('LoPATV_TV')

Result of LoPATV filtering

txt = sprintf('Estimated signal (LoPATV) (RMSE = %.3f)', rmse_lopatv);

figure(2), clf
subplot(2,1,1)
plot(x + p, 'color','black')
title(txt)
printme('LoPATV')

Display with noisy data

txt = sprintf('Estimated signal (LoPATV) (RMSE = %.3f)', rmse_lopatv);

figure(2), clf
subplot(2,1,1)
plot(n, y, '.r', n, x + p, 'black')
title(txt)
printme('LoPATV_color')

Enhanced LoPATV

Lp minimization

mu0 = 50;
mu = .1;
Nit = 200;

lambda = 8;
P = 40;     % P - block overlap
L = 50;     % L - block length

deg = 2;    % deg - degree of polynomial
pow = 0.7;
E = 0.05;

[x, p, cost] = lopatv_Lp(y, L, P, deg, lambda, Nit, mu0, mu, pow, E);

rmse_lopatv_Lp = sqrt(mean((s - x - p).^2));

fprintf('Enhanced LoPATV: RMSE = %.2e\n', rmse_lopatv_Lp)

% Display result
txt = sprintf('Estimated signal (Enhanced LoPATV) (RMSE = %.3f)', rmse_lopatv_Lp);
figure(1)
clf
subplot(2,1,1)
plot(x + p, 'black')
title(txt)
printme('LoPATV_Lp')
Enhanced LoPATV: RMSE = 8.20e-02

Display result with noisy data

figure(1)
clf
subplot(2,1,1)
plot(n, y,'.r', n, x + p, 'black')
title(txt)
printme('LoPATV_Lp_color')

Piecewise constant component (TV component)

figure(2)
clf
subplot(2,1,1)
plot(x, 'black')
title('Calculated TV component (Enhanced LoPATV)');
subplot(2,1,2)
stem(abs(diff(x)),'marker','none', 'color','black')
title('First-order difference');
printme('LoPATV_TV')