resonance_demo: 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 approximation approximation: 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 A1 = 0.35; % lam1, lam2: regularization parameters A2 = 0.40; mu = 0.1; % SALSA parameter % 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) x = x + 0.1 * randn(1,N); 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_noisy_1_%s',str)); % Display signals figure(1), clf subplot(4,1,1) plot(t,x) title('(NOISY) TEST SIGNAL') box off ylim([-A A]) xlim([0 N]) xlabel(xlabel_txt) case 2 % Speech signal x = load('speech1.txt'); fs = 16000; x = x(:)'; N = 2^11; % length(x) t = (0:N-1)/fs; % time axis x = x + 0.02 * randn(1,N); % add noise 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 A1 = 0.1; % lam1, lam2: regularization parameters A2 = 0.1; mu = 0.5; % SALSA parameter % Function for printing a figure to a file printme = @(str) print('-dpdf',sprintf('figures/resonance_demo_noisy_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.318e-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(2), 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')
Peform MCA (decomposition)
now1 = ComputeNow(N,Q1,r1,J1,'radix2'); now2 = ComputeNow(N,Q2,r2,J2,'radix2'); lam1 = A1*now1; lam2 = A2*now2; % [y1,y2,w1s,w2s,costfn] = dualQd(x,Q1,r1,J1,Q2,r2,J2,lam1,lam2,mu,Nit,'plots'); [y1,y2,w1s,w2s,costfn] = dualQd(x,Q1,r1,J1,Q2,r2,J2,lam1,lam2,mu,Nit); % [y1,y2,w1s,w2s] = dualQd(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 = 3.042e-01
Compute cost function
cost = sum(abs(x - y1 - y2).^2); 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); % check consistency between 'cost' and final value of 'costfn' cost_err = costfn(end) - cost; fprintf('Cost function calculation error: %d\n', cost_err)
Objective function: 2.730081e+00 Cost function calculation error: 0
Display cost function versus iteration
The decomposition is based on minimizing a cost function by an iterative optimization algorithm (dualQd). 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(3) clf plot(1:Nit, costfn) xlabel('ITERATION') xlim([0 Nit]) title('COST FUNCTION') grid box off printme('Cost_Function')
Display component signals
The high Q-factor (high-resonance) and low Q-factor (low-resonance) components obtained by minimizing the cost function.
figure(4) clf subplot(4,1,1) plot(t,x) axis([0 N/fs -A A]) title('(NOISY) 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')
Print Information
file_name = sprintf('figures/resonance_demo_noisy_%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,' lambda1 = %.3e \n lambda2 = %.3e \n mu = %.3e\n', lam1, lam2, mu); fprintf(fid,' mu = %.3e\n', mu); fprintf(fid,' iterations = %d\n',Nit); fclose(fid);