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
