demo2: group-sparse thresholding, 2D data (speech spectrogram)

This demo uses 2D overlapping group-sparse (OGS) thresholding for speech denoising. The 2D OGS method is applied to the spectrogram of a noisy speech signal.

This demo illustrates the non-convex form of OGS, with two different non-convex penalties: atan and MCP.

Po-Yu Chen and Ivan Selesnick, August 2013. Revised May 2016.

Reference: P. Y. Chen and I. W. Selesnick, "Group-Sparse Signal Denoising: Non-Convex Regularization, Convex Optimization," IEEE Transactions on Signal Processing, vol. 62, no. 13, pp. 3464-3478, July 1, 2014. doi: 10.1109/TSP.2014.2329274

Set up
Load speech signal
Make noisy signal
Define STFT function
STFT of noisy speech
OGS group-size parameters
Denoising with OGS
OGS with optional arguments
Soft-thresholding

Set up

clear
close all

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

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

Load speech signal

[s, Fs] = audioread('arctic_a0001.wav');

s = s(:)';              % Convert to row vector

N = length(s);          % N : signal length
T = N/Fs;               % T : signal duration (sec)

Make noisy signal

% y : speech signal + noise
sigma = 0.025;
y = s + sigma * randn(1, N);

Define STFT function

Nf = 2^9;                       % FFT length
F = @(y) stft(y, Nf);           % STFT
invF = @(c) istft(c, Nf, N);    % inverse STFT

Y = F(y);    % Y : STFT of noisy speech signal
S = F(s);    % X : STFT of clean speech signal

STFT of noisy speech

dBlim = [-50 0];

figure(1), clf
displaySTFT(Y, Fs, T, dBlim)
title('Noisy data')
print -dpdf figures/demo2_noisy

OGS group-size parameters

% (K1, K2) : Two-dimensional group size in STFT domain
K1 = 8;     % K1 : number of spectral samples
K2 = 2;     % K2 : number of temporal samples

Denoising with OGS

Use the function ogs2 This is non-convex OGS with the default penalty 'atan'.

lam = 0.006;
X = ogs2(Y, K1, K2, lam);    % OGS algorithm
x = real(invF(X));           % inverse STFT

% print RMSE
fprintf('OGS thresholding, atan penalty (default): RMSE = %.4f\n', rmse(s - x))

% Spectrogram
figure(2), clf
displaySTFT(X, Fs, T, dBlim)
title('OGS thresholding, atan penalty (default)');
print -dpdf figures/demo2_ogs

% sound(x, Fs);         % Uncomment to hear processed signal.

OGS thresholding, atan penalty (default): RMSE = 0.0120

OGS with optional arguments

Use ogs2 function with optional arguments to specify penalty and rho parameter (rho between 0 and 1). This is non-convex when rho > 0.

lam = 0.006;
pen = 'mcp';
rho = 0.8;
X = ogs2(Y,  K1, K2, lam, pen, rho);  % OGS algorithm
x = real(invF(X));                    % inverse STFT

% print RMSE
fprintf('OGS thresholding, MCP penalty: RMSE = %.4f\n', rmse(s - x))

% Spectrogram
figure(3), clf
displaySTFT(X, Fs, T, dBlim)
title('OGS thresholding, mcp penalty' )
print -dpdf figures/demo2_ogs_mcp

% sound(x, Fs);         % Uncomment to hear processed signal.

OGS thresholding, MCP penalty: RMSE = 0.0123

Soft-thresholding

Soft-thresholding does not induce group sparsity

lam = 0.05;
X = soft(Y, lam);       % soft-thresholding
x = real(invF(X));      % inverse STFT

% print RMSE
fprintf('Soft thresholding: RMSE = %.4f\n', rmse(s - x))

% Display spectrogram
figure(1), clf
displaySTFT(X, Fs, T, dBlim)
title('Soft thresholding')
print -dpdf figures/demo2_soft

% sound(x, Fs);         % Uncomment to hear processed signal.

Soft thresholding: RMSE = 0.0230