Example2 : Sparse deconvolution of speech

Apply sparse deconvolution to a segment of voiced speech

Contents

Start

clear
close all
printme = @(filename) print('-dpdf', sprintf('figures/Example2_%s', filename));
load mtlb

figure(1)
clf
subplot(2,1,1)
plot(mtlb)

[bb,aa] = butter(2, 0.05, 'high');
z = filtfilt(bb, aa, mtlb);

subplot(2,1,2)
plot(z)
xlim([1000 1600])
box off
y = z(1000:1600);
N = length(y);
n = 0:N-1;
t = n/Fs;

ylim1 = [-3 3];

figure(2)
clf
subplot(2,1,1)
plot(t, y)
box off
ylim(ylim1)
xlim([0 N/Fs])
xlabel('Time (sec)')
title('Speech waveform [y]')

orient landscape
printme('waveform_y')

Find AR model

a = arburg(y, 12);          % AR model from data. Order 12

h = filter(1, a, [1 zeros(1,100)]);         % impulse response of filter

figure(1)
clf
plot(0:100, h)
box off
xlabel('Time (samples)')
title('Impulse response h [AR model]')

printme('impulse_response')

Pole-zero diagram

figure(1)
clf
[zh, ph, lh] = zplane(1, a);
xlabel(' ')
ylabel(' ')
title('Poles/Zeros')
box off
MS = 6;
set(zh, 'markersize', MS)
set(ph, 'markersize', MS)
set(lh(2), 'fontsize', 7)

printme('pole_zero')

Sparse deconvolution

lam = 0.8;                                      % lam : regularization parameter
Nit = 100;                                      % Nit : number of iterations
[x, cost] = deconvL1(y, lam, 1, a, Nit);        % Run sparse deconvolution algorithm

figure(2)
clf
subplot(2,1,1)
plot(1:Nit, cost, '.-')
box off
title('Cost function history');
xlabel('Iteration')
xlim([10 Nit])

Display sparse input signal

figure(2)
clf
subplot(2,1,1)
plot(t, x)
xlim([0 N/Fs])
xlabel('Time (sec)')
box off
title('Sparse deconvolution [x]')
axnote(sprintf('\\lambda = %.2f', lam))

Approximate reconstruction

Approximate reconstruction of original speech waveform

r = filter(1, a, x);

subplot(2,1,2)
plot(t, r)
xlim([0 N/Fs])
xlabel('Time (sec)')
box off
ylim(ylim1)
title('Approximate reconstruction [r = h conv x]')

orient landscape
printme('sparse_deconv')

Exact deconvolution

using A(z)

g = filter(a, 1, y);

% verify:
err = y - filter(1, a, g);
max(abs(err))

figure(1)
clf
subplot(2,1,1)
plot(t, g)
xlim([0 N/Fs])
xlabel('Time (sec)')
box off
title('Exact deconvolution (inverse filter)')

orient landscape
printme('exact_deconv')
ans =

   1.6653e-15