Example: Signal separation using dual basis pursuit (Example 1)

Decompose a signal (spikes plus sinusoids) in to two signals: spikes alone and sinusoids alone. The example illustrated 'morphological component analysis' (MCA) using dual basis pursuit.

Ivan Selesnick
NYU-Poly
selesi@poly.edu
2011

Contents

Start

clear
close all
% I = sqrt(-1);

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

Create signal

The simulated signal contains two spikes and two sinusoids.

M = 150;                                % M : signal length
m = (0:M-1)';
y = cos(0.1*pi*m) + cos(0.15*pi*m);     % sinusoids
y = y/1.5;
y = y + 3*(m == 50);                    % spikes
y = y - 3*(m == 70);
ymax = max(abs(y));

figure(1)
clf
subplot(2,1,1)
plot(m, y);
xlim([0 M])
ylim([-ymax ymax])
box off
title('Signal');
xlabel(' ')

printme('signal')

Define the two Parseval transforms

Define tramsforms A1 and A2 such that A1*A1T = p1*Identity, and A2*A2T = p2*Identity.

% A1 : identity
A1 = @(x) x;
A1T = @(x) x;
p1 = 1;

% A2 : over-sampled DFT (normalized so that Parseval frame constant is 1.)
N = 256;
A2 = @(x) A(x,M,N)/sqrt(N);
A2T = @(x) AT(x,M,N)/sqrt(N);
p2 = 1;

Verify Parseval energy identity

E = sum(abs(y(:)).^2);
c1 = A1T(y);
c2 = A2T(y);
E1 = sum(abs(c1(:)).^2);
E2 = sum(abs(c2(:)).^2);
fprintf('Signal energy = %.3e\n', E)
fprintf('Transform 1 energy = %.3e\n', E1/p1)
fprintf('Transform 2 energy = %.3e\n', E2/p2)

% The energy is the same with both transforms (and the same as the signal energy).
Signal energy = 8.197e+01
Transform 1 energy = 8.197e+01
Transform 2 energy = 8.197e+01

Verify perfect reconstruction property

err = p1*y - A1(c1);
fprintf('Transform A1 : maximum reconstruction error = %g\n', max(abs(err(:))))

err = p2*y - A2(c2);
fprintf('Transform A2 : maximum reconstruction error = %g\n', max(abs(err(:))))
Transform A1 : maximum reconstruction error = 0
Transform A2 : maximum reconstruction error = 4.44089e-16

Peform signal separation using dual BP

Use the command 'dualBP' (dual basis pursuit) to separate the signal into two distinct components.

theta = 0.6;                % theta : trade-off parameter
Nit = 100;                  % Nit : number of iterations
mu1 = 1.5;                  % mu1, mu2 : ADMM parameters
mu2 = 1.0;

[y1,y2,c1,c2,costfn] = dualBP(y, A1, A1T, p1, A2, A2T, p2, theta, 1-theta, mu1, mu2, Nit);

Display cost function

figure(2)
clf
it1 = 1;
plot(1:Nit, costfn)
xlim([0 Nit])
box off
title('Cost function')
xlabel('Iteration')

Calculate signal components

The two components y1 and y2 can be found by applying the transforms to the coefficients c1 and c2 produced by the dual basis pursuit algorithhm.

y1 = A1(c1);
y2 = A2(c2);

% Verify that y = y1 + y2
fprintf('Maximum of residual = %g\n', max(abs(y - y1 - y2)))
Maximum of residual = 4.44089e-16

Display signal components obtained using dual BP

figure(3)
clf
subplot(2,1,1)
plot(m, real(y1))
xlim([0 M])
title('Component 1');
box off

subplot(2,1,2)
plot(m, real(y2))
xlim([0 M])
box off
title('Component 2');
xlabel(' ')

printme('components')

It can be seen that the coefficients are sparse.

figure(4)
clf
subplot(2,1,1)
stem(c1, 'marker', 'none')
box off
title('Coefficients c1');

subplot(2,1,2)
stem(abs(c2), 'marker', 'none')
title('Coefficients c2');
box off
xlim([0 N])
xlabel(' ')
printme('coefficients')