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