GMC_demo_1

Demonstration of the generalized MC (GMC) penalty for sparse-regularized linear least squares. The GMC is designed to maintain the convexity of the cost function to be minimized. In this example, sparse regularization is used to estimate a signal (with frequency-domain sparsity) form a noisy observation.

Demonstration software for the paper: I. Selesnick, 'Sparsity amplified', IEEE Int. Conf. Acoust., Speech, Signal Processing (ICASSP), Mar. 2017.

Ivan Selesnick
May 2016
Revised: July 2016

Contents

Start

clear

% Set random number generator for reproducibility
rng('default')
rng(1)

rmse = @(err) sqrt(mean(abs(err(:)).^2));

Define transform

The transform A is wide-matrix inverse DFT. A is of size 100 x 256. A satisfies A A' = I

M = 100;    % M: signal length
N = 256;    % N: transform length

truncate = @(x, M) x(1:M);

% A: inverse DFT (normalized for Parsevals property)
AH = @(x) fft(x, N)/sqrt(N);
A = @(X) truncate(ifft(X), M) * sqrt(N);
normA = sqrt(M/N);

Verify properties of A

% Verify that A AH = I
y = randn(M, 1);
err = A(AH(y)) - y;
max(abs(err))

% Verify further:
Amtx = dftmtx(N)'/sqrt(N); Amtx = Amtx(1:M, :);

r = randn(M, 1);
err = Amtx' * r - AH(r);
max(abs(err))
v = rand(size(AH(r)));
err = Amtx * v - A(v);
max(abs(err))
% max(eig(Amtx'*Amtx))

ans =
   4.4409e-16
ans =
   6.6613e-16
ans =
   7.1054e-15

Create data

Sinusoids plus white noise

f1 = 0.1;
f2 = 0.22;
m = (0:M-1)';
s = 2*cos(2*pi*f1*m) + sin(2*pi*f2*m);    % s : true signal

sigma = 1;           % sigma : noise standard deviation
w = randn(M, 1);     % w : white Gaussian noise
y = s + sigma * w;   % y : signal plus noise

figure(1)
clf
subplot(2, 1, 1)
plot(m, s)
title('Noise-free data')
xlabel('Time')

subplot(2, 1, 2)
plot(m, y)
title('Noisy data')
xlabel('Time')

L1 norm

Least squares with L1 norm regularization

% Algorithm parameters
beta_L1 = 1.8;
lam_L1 = beta_L1 * sigma * normA;

c_L1 = srls_L1(y, A, AH, 1, lam_L1);
x_L1 = A(c_L1);

rmse_L1 = rmse(x_L1 - s)

rmse_L1 =
    0.4435

GMC penalty

Least squares with GMC regularization

beta_GMC = 2.6;
lam_GMC = beta_GMC * sigma * normA;
gamma = 0.8;

c_GMC = srls_GMC(y, A, AH, 1, lam_GMC, gamma);
x_GMC = A(c_GMC);

rmse_GMC = rmse(x_GMC - s)

rmse_GMC =
    0.2221

Plot sparse coefficients

figure(1)
clf
plot(1:N, abs(c_GMC), 'o', 1:N, abs(c_L1), 'x')
legend('GMC', 'L1', 'location', 'north')
title('Sparse Fourier coefficients')
xlim([0 N])

Plot solutions

f = (0:N-1)/N;

figure(2)
clf

% Noise-free signal
subplot(3, 2, 1)
plot(m, real(s))
title('Noise-free signal');
ylim1 = [-1 1]*max(abs(y)) * 1.1;
ylim(ylim1)
xlabel('Time (n)')

% Noisy data
subplot(3, 2, 2)
plot(m, real(y))
title(sprintf('Noisy signal, \\sigma = %.2f', sigma));
ylim(ylim1)
xlabel('Time (n)')

subplot(3, 2, 5)
plot(f, abs(AH(y)))
title('FFT of noisy signal');
xlabel('Frequency (normalized)')
xlim([0 0.5])
ylim([0 8])

% L1-norm
subplot(3, 2, 3)
plot(m, x_L1)
title(sprintf('Denoising using L1 norm, RMSE = %.3f', rmse(x_L1 - s)))
ylim(ylim1)
xlabel('Time (n)')

% GMC-penalty
subplot(3, 2, 4)
plot(m, x_GMC)
title(sprintf('Denoising using GMC penalty, RMSE = %.3f', rmse(x_GMC - s)))
ylim(ylim1)
xlabel('Time (n)')

% Sparse coefficients
subplot(3, 2, 6)
k_L1 = abs(c_L1) > eps;
k_GMC = abs(c_GMC) > eps;
ph_L1 = plot(f(k_L1), abs(c_L1(k_L1)), 'x');
hold on
ph_GMC = plot(f(k_GMC), abs(c_GMC(k_GMC)), 'o', 'markersize', 4);
legend([ph_GMC, ph_L1], 'GMC', 'L1 norm')
title('Optimized coefficients');
xlabel('Frequency (normalized)')
ylim([0 15])
xlim([0 0.5])

orient landscape
print -dpdf figures/GMC_demo

Published with MATLAB® 7.12