Example: Basis pursuit denoising (BPD) with the STFT

Estimation of pulse in white Gaussian noise. The STFT is used.

Contents

Misc

clear
close all

I = sqrt(-1);

dB = @(x) 20 * log10(abs(x));

Make signal

N = 100;                   % N : signal length
n = 0:N-1;

f = 0.15;
K = 20;                    % K : duration of pulse
M = 55;                    % M : start-time of pulse
A = 3;
s = zeros(size(n));
s(M+(1:K)) = A*sin(2*pi*f*(1:K)) .* hamming(K)';
% s : pulse

figure(1)
clf
plot(n, s)
title('Noise-free signal [s]')

Define transform

R = 16;         % R : frame length
Nfft = 32;      % Nfft : FFT length in STFT (Nfft >= R)

[A, AH, normA] = MakeTransforms('STFT', N, [R Nfft]);
Reconstruction error = 0.000000
Energy ratio = 1.000000

Validate normA calculation

% imp : an impulse in transform domain
imp = AH(zeros(1, N));
n0 = 10;
k0 = 5;
imp(n0, k0) = 1;

[normA norm(A(imp))]   % should be same
ans =

    0.5000    0.5000

Make noisy data

randn('state', 0); % set state for reproducibility of example

sigma = 0.5;
y = s + sigma * randn(size(s));

Display noisy data and its STFT

figure(1)
clf

subplot(2, 1, 1)
plot(n, y)
title('Noisy data [y]')


AHy = AH(y);

tt = R/2 * ( (0:size(AHy, 2)) - 1 );    % tt : time axis for STFT

subplot(2, 1, 2)
imagesc(tt, [0 0.5], dB(AHy(1:Nfft/2+1, :)), max(dB(AHy(:))) + [-50 0])
axis xy
xlim([0 N])
ylim([0 0.5])
title('STFT of Noisy data')
ylabel('Frequency')
xlabel('Time')

cm = colormap(gray);                    % cm : colormap for STFT
cm = cm(end:-1:1, :);
colormap(cm)

Solve BPD

beta = 2.5;
lam = beta * sigma * normA;

mu = 0.1;
Nit = 100;

[c, cost] = BPD(y, A, AH, lam, mu, Nit);

x = A(c);

% Display cost function history to observe convergence of algorithm.
figure(1)
clf
plot(cost)
title('Cost function history')
xlabel('Iteration')

Display denoised signal

cmax = max(abs(c(:)));

figure(1)
clf

subplot(2, 1, 1)
plot(n, real(x))
title('Denoised signal')
ylabel('real( Ac )')
ch = colorbar;
set(ch, 'visible', 'off')

subplot(2, 1, 2)
imagesc(tt, [0 0.5], dB(c(1:Nfft/2+1, :)), dB(cmax) + [-50 0])
% imagesc(tt, [0 1], dB(c), dB(cmax) + [-50 0])
axis xy
xlim([0 N])
ylim([0 0.5])
% title('BPD STFT coefficients');
title(sprintf('BPD STFT coefficients (beta = %.2f)', beta));
ylabel('Frequency')
xlabel('Time')
colormap(cm)
ch = colorbar;

Optimality scatter plot

I = sqrt(-1);
g = (1/lam) * AH(y - A(c)) .* exp(-I*angle(c));

figure(2)
clf
plot3( real(g(:)), imag(g(:)), abs(c(:)), '.')
zlabel('abs( c )')
xlabel('Re(A^H(y - A c) e^{-j\anglec})/\lambda')
ylabel('Im(A^H(y - A c) e^{-j\anglec})/\lambda')
title('Optimality scatter plot')
grid off
theta = linspace(0, 2*pi, 200);
line( cos(theta), sin(theta), 'color', [1 1 1]*0.5)
line([1 1], [0 0], [0 cmax],  'color', [1 1 1]*0.5)
xm = 1.2;
xlim([-1 1]*xm)
ylim([-1 1]*xm)
line([-1 1]*xm, [1 1]*xm, 'color', 'black')
line([1 1]*xm, [-1 1]*xm, 'color', 'black')

% Animate: vary view orientation
az = -36;
for el = [40:-1:5]
    view([az el])
    drawnow
end
for az = [-36:-3]
    view([az el])
    drawnow
end

Save figure to file

MyGraphPrefs('on')

figure(1)
clf

subplot(4, 1, 1)
line(n, y)
mytitle('Noisy data');
ylabel('y')
xlabel('Time')
ylim([-3 3])

subplot(4, 1, 2)
imagesc(tt, [0 0.5], dB(AHy(1:Nfft/2+1, :)), max(dB(AHy(:))) + [-50 0])
% imagesc(tt, [0 1], dB(abs(AH(y))), dB(cmax) + [-50 0])
axis xy
colormap(cm)
xlim([0 N])
ylim([0 0.5])
mytitle('STFT of noisy data');
xlabel('Time')
ylabel('Frequency')
box off

subplot(4, 1, 3)
imagesc(tt, [0 0.5], dB(c(1:Nfft/2+1, :)), dB(cmax) + [-50 0])
% imagesc(tt, [0 1], dB(c), dB(cmax) + [-50 0])
axis xy
colormap(cm)
xlim([0 N])
ylim([0 0.5])
mytitle(sprintf('BPD STFT coefficients (beta = %.2f)', beta));
xlabel('Time')
ylabel('Frequency')
box off


subplot(4, 1, 4)
line(n, real(A(c)))
mytitle('Denoised signal');
ylabel('real( Ac )')
xlabel('Time')
ylim([-3 3])


% print figure to pdf file
orient tall
print -dpdf figures/BPD_example_stft

% print figure to eps file
set(gcf, 'PaperPosition', [1 1 4 7.5])
print -deps figures_eps/BPD_example_stft

MyGraphPrefs('off')

Animate: vary lambda

for beta = 0.1:0.1:3

    lam = beta * sigma * normA;
    [c, cost] = BPD(y, A, AH, lam, mu, Nit);
    x = A(c);

    figure(10)
    clf

    subplot(2, 1, 1)
    imagesc(tt, [0 0.5], dB(c(1:Nfft/2+1, :)), dB(cmax) + [-50 0])
%     imagesc(tt, [0 1], dB(abs(c)), [-20 20])
    axis xy
    xlim([0 N])
    ylim([0 0.5])
    title(sprintf('BPD STFT coefficients : \\lambda = \\beta\\sigma ||a||, \\beta = %.2f', beta));
    ylabel('Frequency')
    xlabel('Time')
    colormap(cm)

    subplot(2, 1, 2)
    plot(n, real(x))
    ylim([-1 1]*4)
    title(sprintf('Denoised signal, \\beta = %.2f', beta));
    ylabel('real( Ac )')
    xlabel('Time')

    drawnow

end