MATLAB Examples

resonance_demo_eq: Resonance-based signal decomposition

Example illustrating the decomposition of a test signal into a high-resonance component and a low-resonance component (into an oscillatory component and a transient component).

This version is based on exact equality x = x1 + x2

Reference: Resonance-Based Signal Decomposition: A New Sparsity-Enabled Signal Analysis Method. Signal Processing, 2010, doi:10.1016/j.sigpro.2010.10.018 I. W. Selesnick.

Contents

Miscellaneous

close all
clear

Set example parameters

% Uncomment one of the following two lines
% Example_Num = 1;              % Artificial signal
Example_Num = 2;                % Speech waveform

switch Example_Num

    case 1

        % Artifical signal

        % High Q-factor wavelet transform parameters
        Q1 = 3;
        r1 = 3;
        J1 = 19;

        % Low Q-factor wavelet transform parameters
        Q2 = 1;
        r2 = 3;
        J2 = 9;

        % Set MCA parameters
        Nit = 100;                  % Number of iterations
        mu = 0.5;                   % SALSA parameter
        theta = 0.47;

        % Make test signal
        N = 2^9; % N = 512;         % length of signal x is power of 2
        x1 = test_signal(5);        % High-resonance signal
        x2 = test_signal(1);        % Low-resonance signal
        x = x1 + x2;                % Test signal (sum of high and low resonances)
        t = (0:N-1);                % time axis
        fs = 1;
        xlabel_txt = 'TIME (SAMPLES)';
        A = 2;                      % y-axis extent for plots

        % Function for printing a figure to a file
        printme = @(str) print('-dpdf',sprintf('figures/resonance_demo_1_%s',str));

        % Display signals
        figure(1), clf
        subplot(4,1,1)
        plot(t,x)
        title('TEST SIGNAL')
        box off
        ylim([-A A])
        xlim([0 N])

        subplot(4,1,2)
        plot(t,x1)
        box off
        axis([0 N -A A])
        title('HIGH-RESONANCE COMPONENT')

        subplot(4,1,3)
        plot(t,x2)
        box off
        axis([0 N -A A])
        title('LOW-RESONANCE COMPONENT')
        xlabel(xlabel_txt)

        orient tall
        printme('Test_Signals')

    case 2

        % Speech signal
        x = load('speech1.txt');
        fs = 16000;
        x = x(:)';
        N = 2^11;               % length(x)
        t = (0:N-1)/fs;         % time axis

        xlabel_txt = 'TIME (SECONDS)';
        A = 0.3;

        % High Q-factor wavelet transform parameters
        Q1 = 4;
        r1 = 3;
        J1 = 30;

        % Low Q-factor wavelet transform parameters
        Q2 = 1;
        r2 = 3;
        J2 = 10;

        % Set MCA parameters
        Nit = 100;          % Number of iterations
        mu = 5.0;           % SALSA parameter
        theta = 0.5;

        % Function for printing a figure to a file
        printme = @(str) print('-dpdf',sprintf('figures/resonance_demo_2_%s',str));

        % Display signal
        figure(1), clf
        subplot(4,1,1)
        plot(t,x)
        title('SEECH WAVEFORM')
        box off
        ylim([-A A])
        xlim([0 N/fs])
        xlabel(xlabel_txt)

        orient landscape
        printme('Speech_Waveform')

end

Verify perfect reconstruction

Check that transforms with selected parameters satisfy perfect reconstruction. (It is not really needed; just to double check.)

w1 = tqwt_radix2(x, Q1, r1, J1);
y = itqwt_radix2(w1, Q1, r1, N);

fprintf('Transform 1 reconstruction error = %4.3e\n', max(abs(y-x)))

w2 = tqwt_radix2(x, Q2, r2, J2);
y = itqwt_radix2(w2, Q2, r2, N);

fprintf('Transform 2 reconstruction error = %4.3e\n', max(abs(y-x)))
Transform 1 reconstruction error = 1.110e-16
Transform 2 reconstruction error = 1.388e-16

Plot wavelets at final levels

% Compute wavelets
wlets1 = ComputeWavelets(N,Q1,r1,J1,'radix2');
wlets2 = ComputeWavelets(N,Q2,r2,J2,'radix2');

figure(4), clf
subplot(2,1,1)
plot(t, wlets1{J1})
title(sprintf('HIGH Q-FACTOR WAVELET AT LEVEL %d',J1))
xlim([0 N/fs])
ylim(1.2*(max(abs(wlets1{J1})))*[-1 1])
box off
xlabel(xlabel_txt)
subplot(2,1,2)
plot(t, wlets2{J2})
title(sprintf('LOW Q-FACTOR WAVELET AT LEVEL %d',J2))
xlim([0 N/fs])
ylim(1.2*(max(abs(wlets2{J2})))*[-1 1])
box off
xlabel(xlabel_txt)

printme('Wavelets')

Display energy distribution

It can be useful to know how the energy of a signal is distributed accross the subbands. We display the energy in each subband. Because the transform has the Parseval property the distribution of the energy accross the subbands reflects the frequency content of the signal.

figure(5)
clf
subplot(2,1,1)
PlotEnergy(w1);
title('DISTRIBUTION OF SIGNAL ENERGY: HIGH Q-FACTOR TRANSFORM')
subplot(2,1,2)
PlotEnergy(w2);
title('DISTRIBUTION OF SIGNAL ENERGY: LOW Q-FACTOR TRANSFORM')
printme('Energy')

Peform MCA (decomposition)

now1 = ComputeNow(N,Q1,r1,J1,'radix2');
now2 = ComputeNow(N,Q2,r2,J2,'radix2');
lam1 = theta * now1;
lam2 = (1-theta) * now2;

% [y1,y2,w1s,w2s,costfn] = dualQ(x,Q1,r1,J1,Q2,r2,J2,lam1,lam2,mu,Nit,'plots');
[y1,y2,w1s,w2s,costfn] = dualQ(x,Q1,r1,J1,Q2,r2,J2,lam1,lam2,mu,Nit);

res = x - y1 - y2;  % residual

rel_res = sqrt(mean(res.^2))/sqrt(mean(x.^2));

fprintf('relative RMS reconstruction error = %4.3e\n',rel_res);

% check
% max(abs(y1-itqwt_radix2(w1s,Q1,r1,N)))        % should be zero
% max(abs(y2-itqwt_radix2(w2s,Q2,r2,N)))        % should be zero
relative RMS reconstruction error = 4.354e-16

Compute cost function

cost = 0;
for j = 1:J1+1
    cost = cost + lam1(j)*sum(abs(w1s{j}));
end
for j = 1:J2+1
    cost = cost + lam2(j)*sum(abs(w2s{j}));
end

fprintf('Objective function: %e\n', cost);
Objective function: 1.259349e+01

Display cost function versus iteration

The decomposition is based on minimizing a cost function by an iterative optimization algorithm (dualQ). It can be useful to inspect the value of the cost function versus the iteration number to see if the algorithm has converged satisfactorily.

figure(6)
clf
plot(1:Nit, costfn)
xlabel('ITERATION')
xlim([0 Nit])
title('COST FUNCTION')
grid
box off
printme('Cost_Function')

% check consistency between 'cost' and final value of 'costfn'
cost_err = costfn(end) - cost;
fprintf('Cost function calculation error: %d\n', cost_err)
Cost function calculation error: 0

Display component signals

The high Q-factor (high-resonance) and low Q-factor (low-resonance) components obtained by minimizing the cost function.

figure(7)
clf

subplot(4,1,1)
plot(t,x)
axis([0 N/fs -A A])
title('SIGNAL')
box off

subplot(4,1,2)
plot(t,y1)
axis([0 N/fs -A A])
title('HIGH-RESONANCE COMPONENT')
box off

subplot(4,1,3)
plot(t,y2)
axis([0 N/fs -A A])
title('LOW-RESONANCE COMPONENT')
box off

subplot(4,1,4)
plot(t, res)
axis([0 N/fs -A A])
title('RESIDUAL')
box off

xlabel(xlabel_txt)

orient tall
if Example_Num == 2
    orient landscape
end
printme('Components')

Display energy distribution

Display the distribution of energy of the high Q-factor and low Q-factor representations obtained by minimizing the cost function.

figure(8)
clf
subplot(2,1,1)
PlotEnergy(w1s);
title('DISTRIBUTION OF SIGNAL ENERGY: HIGH Q-FACTOR TRANSFORM')
subplot(2,1,2)
PlotEnergy(w2s);
title('DISTRIBUTION OF SIGNAL ENERGY: LOW Q-FACTOR TRANSFORM')
printme('Energy_Sparse')

Print Information

file_name = sprintf('figures/resonance_demo_%d_info.txt',Example_Num);
fid = fopen(file_name,'w');

fprintf(fid,'TRANSFORM TYPE: tqwt_radix2\n');
fprintf(fid,' TQWT 1: Q = %4.2f, r = %4.2f, levels = %d\n',Q1,r1,J1);
fprintf(fid,' TQWT 2: Q = %4.2f, r = %4.2f, levels = %d\n\n',Q2,r2,J2);
fprintf(fid,'SALSA PARAMETERS: \n');
fprintf(fid,' mu = %.3e\n', mu);
fprintf(fid,' iterations = %d\n',Nit);

fclose(fid);