Sparse representation with the TQWT
Example illustrating sparse signal representation/approximation using the tunable Q-factor wavelet transform (TQWT). The first part: sparse signal representation (basis pursuit). The second part: sparse signal approximation (basis pursuit densoising).
Reference: 'Wavelet Transform with Tunable Q-Factor' http://taco.poly.edu/selesi/TQWT/ Ivan Selesnick, selesi@poly.edu Polytechnic Institute of NYU November 2010
Contents
- Miscellaneous
- Set example parameters
- Display signal
- Plot wavelets for several scales
- Plot wavelet at final level and low-pass scaling function
- Plot subbands
- Compute energy in each subband
- Sparse wavelet representation (Basis Pursuit)
- Compute cost function
- Display cost function versus iteration
- Plot sparse subbands
- Reconstruction error
- Compute energy in each subband
- Sparse wavelet approximation (Basis Pursuit Denoising)
- Compute cost function
- Display cost function versus
- Display test signal before and after denoising
- Plot sparse subbands
- Residual
- Write information to file
Miscellaneous
clear
close all
Set example parameters
% Uncomment one of the following lines Example_Num = 1; % Artificial signal % Example_Num = 2; % Speech waveform switch Example_Num case 1 % Set wavelet parameters Q = 3.5; r = 3.0; L = 20; % number of levels L1 = 6; x = test_signal(2); % Make test signal N = length(x); t = (0:N-1); % time axis fs = 1; % sampling frequency xlabel_txt = 'TIME (SAMPLES)'; A = 2; printme = @(str) print('-dpdf',sprintf('figures/sparsity_demo1_%s',str)); case 2 % Set wavelet parameters Q = 3; r = 3; L = 23; L1 = 10; % Load speech signal x = load('speech2.txt'); fs = 16000; x = x(:)'; N = 2^11; % length(x) t = (0:N-1)/fs; % time axis xlabel_txt = 'TIME (SECONDS)'; A = 0.5; printme = @(str) print('-dpdf',sprintf('figures/sparsity_demo2_%s',str)); end beta = 2/(Q+1); alpha = 1-beta/r;
Display signal
figure(1), clf subplot(2,1,1) plot(t,x) title('TEST SIGNAL') box off ylim([-A A]) xlim([t(1) t(end)]) xlabel(xlabel_txt) % Verify perfect reconstruction w = tqwt_radix2(x, Q, r, L); y = itqwt_radix2(w, Q, r, N); fprintf('Reconstruction error = %4.3e\n', max(abs(y-x)))
Reconstruction error = 5.734e-16
Plot wavelets for several scales
Verify that the wavelets look as expected
figure(2), clf N1 = 200; % Length of signal PlotWavelets(N1,Q,r,4,12); xlabel('TIME (SAMPLES)') orient tall printme('wavelets')
Plot wavelet at final level and low-pass scaling function
wlets = ComputeWavelets(N,Q,r,L,'radix2'); % Compute wavelets figure(3), clf subplot(2,1,1) plot(t, wlets{L}) title(sprintf('WAVELET AT LEVEL %d',L)) xlim([0 N/fs]) ylim(1.2*(max(abs(wlets{L})))*[-1 1]) box off subplot(2,1,2) plot(t, wlets{L+1}) title(sprintf('LOW-PASS SCALING FUNCTION AT LEVEL %d',L)) xlim([0 N/fs]) box off xlabel(xlabel_txt)
Plot subbands
figure(4), clf % PlotSubbands(x,w,Q,r,L1,L,fs,[],'stem'); PlotSubbands(x,w,Q,r,L1,L,fs,'E','stem'); title('SUBBANDS') orient tall printme('subbands')
Compute energy in each subband
It can be useful to know how the energy of a signal is distributed across the subbands. We compute the energy in each subband and display using a bar graph. Because the transform has the Parseval property the distribution of the energy across the subbands reflects the frequency content of the signal.
figure(5)
clf
subplot(3,1,1:2)
e = PlotEnergy(w);
printme('energy')
Sparse wavelet representation (Basis Pursuit)
Sparse signal representation with perfect reconstruction (Basis Pursuit). Use a variant of SALSA to minimize l1-norm of wavelet coefficients providing exact reconstruction.
now = ComputeNow(N,Q,r,L,'radix2'); lambda = now; % Regularization parameters mu = 2.0; % SALSA parameter Nit = 100; % Number of iterations [w2, costfn] = tqwt_bp(x, Q, r, L, lambda, mu, Nit); y = itqwt_radix2(w2, Q, r, N); err = x - y; rel_err = sqrt(mean(abs(err).^2))/sqrt(mean(abs(x).^2)); fprintf('Basis pursuit (BP) relative RMS reconstruction error = %4.3e\n',rel_err);
Basis pursuit (BP) relative RMS reconstruction error = 3.120e-16
Compute cost function
cost = 0; for j = 1:L+1 cost = cost + lambda(j)*sum(abs(w2{j})); end fprintf('BP objective function: %e\n', cost);
BP objective function: 2.221250e+01
Display cost function versus iteration
figure(6), clf it1 = 10; plot(it1:Nit, costfn(it1:Nit)); xlim([0 Nit]) title('SALSA COST FUNCTION (BASIS PURSUIT)') xlabel('ITERATION') box off printme('cost_function_bp') % check consistency between 'cost' and final value of 'costfn' fprintf('BP costfn(end) = %d\n', costfn(end))
BP costfn(end) = 2.221250e+01
Plot sparse subbands
figure(7), clf PlotSubbands(y,w2,Q,r,L1,L,fs,'E','stem'); title('SPARSE SUBBANDS (BASIS PURSUIT)') orient tall printme('subbands_sparse_bp')
Reconstruction error
The reconstruction error is zero, which verifies that the sparse wavelet coefficients are a valid representation of the signal.
figure(8), clf subplot(2,1,1) plot(t, x) xlim([t(1) t(end)]) ylim([-A A]) box off title('TEST SIGNAL') subplot(2,1,2) plot(t, x-y) xlim([t(1) t(end)]) box off title('RECONSTRUCTION ERROR (BASIS PURSUIT)') xlabel(xlabel_txt) printme('recon_error_bp')
Compute energy in each subband
figure(9) clf subplot(3,1,1:2) e2 = PlotEnergy(w2); title('DISTRIBUTION OF SIGNAL ENERGY (BASIS PURSUIT)') printme('energy_bp')
Sparse wavelet approximation (Basis Pursuit Denoising)
Sparse signal representation with l1-norm regularization (Basis Pursuit Denoising). Use SALSA to minimize function: sum((x-invTQWT(w)).^2) + sum(abs((lambda.*w))). This is useful when the signal is noisy.
lambda = 0.5*now; % Regularizaton parameter mu = 0.10; % SALSA parameter Nit = 100; % Number of iterations x2 = x + 0.1*randn(1,N); % Noisy signal [wy, costfn] = tqwt_bpd(x2, Q, r, L, lambda, mu, Nit); y = itqwt_radix2(wy, Q, r, N);
Compute cost function
cost = sum(abs(x2 - y).^2); for j = 1:L+1 cost = cost + lambda(j)*sum(abs(wy{j})); end fprintf('BPD objective function: %e\n', cost);
BPD objective function: 1.507783e+01
Display cost function versus
figure(11), clf it1 = 10; plot(it1:Nit, costfn(it1:Nit)); xlim([0 Nit]) title('SALSA COST FUNCTION (BASIS PURSUIT DENOISING)') xlabel('ITERATION') box off printme('cost_function_bpd') % check consistency between 'cost' and final value of 'costfn' fprintf('BPD costfn(end) = %d\n', costfn(end))
BPD costfn(end) = 1.507783e+01
Display test signal before and after denoising
figure(20), clf subplot(2,1,1) plot(t, x2) xlim([t(1) t(end)]) ylim([-A A]) box off title('NOISY TEST SIGNAL') xlabel(xlabel_txt) subplot(2,1,2) plot(t, y) xlim([t(1) t(end)]) ylim([-A A]) box off title('AFTER BASIS PURSUIT DENOISING') xlabel(xlabel_txt) printme('signals_bpd')
Plot sparse subbands
figure(10), clf PlotSubbands(y,wy,Q,r,L1,L,fs,'E','stem'); title('SPARSE SUBBANDS (BASIS PURSUIT DENOISING)') orient tall printme('subbands_sparse_bpd')
Residual
The residual is an estimate of the noise - it should look like pure noise
figure(11), clf subplot(2,1,1) plot(t, x2-y) xlim([t(1) t(end)]) box off % ylim([-A A]*0.1) title('RESIDUAL (BASIS PURSUIT DENOISING)') xlabel(xlabel_txt) printme('residual_bpd')
Write information to file
file_name = sprintf('figures/sparsity_demo%d_info.txt',Example_Num); fid = fopen(file_name,'w'); fprintf(fid,'TRANSFORM PARAMETERS:\n'); fprintf(fid,'\t Transform: tqwt_radix2.m\n'); fprintf(fid,'\t Q = %4.2f\n\t r = %4.2f\n\t levels = %d\n\n',Q,r,L); fprintf(fid,'SALSA PARAMETERS: \n'); fprintf(fid,'\t mu = %.2e\n\t iterations = %d\n', mu, Nit); fclose(fid);