demo_FFT

Demonstration of the generalized MC (GMC) penalty for sparse-regularized linear least squares. The GMC penalty 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) from a noisy observation.

Demonstration software for the paper: I. Selesnick, 'Sparse Regularization via Convex Analysis' IEEE Transactions on Signal Processing, 2017.

Ivan Selesnick
May 2016
Revised: July 2016, March 2021

Contents

Start

clear
set_plot_defaults

% 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 transform properties

Verify properties of A

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

% Verify A as a matrix
Amtx = dftmtx(N)'/sqrt(N);
Amtx = Amtx(1:M, :);

x = randn(M, 1);
err = Amtx' * x - AH(x);
max(abs(err))
ans =

   4.4409e-16


ans =

   4.5776e-16

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

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

rmse_L1 = rmse(x_L1 - s)

figure(1)
clf
% Signals
subplot(2, 1, 1)
plot(m, real(s), m, x_L1)
title(sprintf('Denoising using L1 norm, RMSE = %.3f', rmse(x_L1 - s)))
legend('True signal', 'L1 norm');
ylim1 = [-1 1]*max(abs(y)) * 1.1;
ylim(ylim1)
xlabel('Time (n)')

% Sparse coefficients
subplot(2, 1, 2)
plot(1:N, abs(c_L1), 'o')
title('Sparse Fourier coefficients (L1 norm)')
xlim([0 N])
xlabel('Frequency (index)')
ylim([0 15])
drawnow

print -dpdf -bestfit figures/demo_FFT_L1norm
rmse_L1 =

    0.4406

GMC penalty

Least squares with GMC regularization

lam_GMC = 2.0;
gamma = 0.8;
c_GMC_ = sparseLS_GMC(y, A, AH, 1, lam_GMC, gamma);

Altnerate algorithm

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

figure(2)
clf
% plot(1:N, abs(c_GMC), 1:N, abs(c_GMC_ver2))
plot( 1:N, abs(c_GMC_), 1:N, abs(c_GMC_ - c_GMC), 'o')
legend('coefficients', 'erorr')

x_GMC = A(c_GMC);

rmse_GMC = rmse(x_GMC - s)
rmse_GMC =

    0.2993

figure(1)
clf
% Signals
subplot(2, 1, 1)
plot(m, real(s), m, x_GMC)
title(sprintf('Denoising using GMC, RMSE = %.3f', rmse(x_GMC - s)))
legend('True signal', 'GMC penalty');
ylim(ylim1)
xlabel('Time (n)')

% Sparse coefficients
subplot(2, 1, 2)
plot(1:N, abs(c_GMC), 'o')
title('Sparse Fourier coefficients (GMC penalty)')
xlim([0 N])
xlabel('Frequency (index)')
ylim([0 15])
drawnow

print -dpdf -bestfit figures/demo_FFT_GMC

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(2, 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(2, 2, 2)
plot(m, real(y))
title(sprintf('Noisy signal, \\sigma = %.2f', sigma));
ylim(ylim1)
xlabel('Time (n)')

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

% Sparse coefficients
subplot(2, 2, 4)
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')
ylim([0 15])
% xlim([0 0.5])
drawnow

orient landscape
print -dpdf -fillpage figures/demo_FFT