MATLAB Examples

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

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);