diff --git a/demos/demo_dequantization.m b/demos/demo_dequantization.m index 06188df..5fc0681 100644 --- a/demos/demo_dequantization.m +++ b/demos/demo_dequantization.m @@ -1,7 +1,46 @@ -% This demo shows how a quantized signal, sparse in the DCT domain, can be dequantized -% solving a convex problem using Douglas-Rachford algorithm +%DEMO_DEQUANTIZATION +% This demo shows how a quantized signal, sparse in the DCT domain, can be dequantized +% solving a convex problem using Douglas-Rachford algorithm +% +% Suppose signal y has been quantized. In this demo we use quantization levels +% that are uniformly spread between the min. and max. value of the +% signal. The resulting signal is y_Q. +% +% The problem can be expressed as +% +% .. argmin_x || x ||_1 s.t. ||Dx - y_Q||_infty <= alpha/2 +% +% .. math:: arg \min_x \| x \|_{1} \text{ s.t. } \|Dx - y_Q\|_\infty \leq \frac{\alpha}{2} +% +% where D is the synthesis dictionary (DCT in our case) and $\alpha$ is the +% distance between quantization levels. The constraint basically +% represents the fact that the reconstructed signal samples must stay +% within the corresponding quantization stripes. +% +% After sparse coordinates are found, the dequantized signal +% is obtained simply by synthesis with the dictionary. +% +% The program is solved using Douglas-Rachford algorithm. We set +% +% * $f_1(x)=||x||_{1}$. Its respective prox is the soft thresholding operator. +% +% * $f_2(x)=i_C$ is the indicator function of the set C, defined as +% .. C = { x | ||Dx - y_Q||_infty <= alpha/2 } +% +% .. math:: C = \{ x | \|Dx - y_Q\|_\infty <= \frac{\alpha}{2} \} +% Its prox is the orthogonal projection onto that set, which is realized +% by entry-wise 1D projections onto the quantization stripes. This is +% realized for all the entries at once by function proj_box. +% +% As an alternative, setting algorithm = 'LP' switches to computing the +% result via linear programming (requires Matlab optimization toolbox). +% +% References: combettes2007douglas + +% Authors: Pavel Rajmic, Pavel Záviška +% Date: August 2015 + -% 2015 Pavel Rajmic, Pavel Záviška clc clear @@ -43,8 +82,8 @@ y = A(x); %or just load a fixed dataset -clear -load('dequantization_dataset_01') +% clear +% load('dequantization_dataset_01') % show coefficients fig_coef = figure; @@ -61,7 +100,7 @@ max_y = max(max(y)); range = max_y - min_y; -quant_step = range / (d-1); %quantization step +quant_step = range / (d-1); %quantization step dec_bounds = (min_y+quant_step/2) : quant_step : (max_y-quant_step/2); %decision boundaries lie in the middle quant_levels = min_y : quant_step : max_y; %quantization levels @@ -172,54 +211,56 @@ info sol = f2.prox(sol,[]); %final projection into the constraints - pause - - %%%%%%%% manual version %%%%%%%%%%%%% - %starting point - DR_y = A(y_quant); - DR_x_old = DR_y; - - relat_change_coefs = 1; - relat_change_obj = 1; - cnt = 1; %iteration counter - obj_eval = []; - -% while relat_change_coefs > 0.00001 - while relat_change_obj > paramsolver.tol - % DR: gamma = 1 - DR_x = f2.prox(DR_y,[]); - obj_eval = [obj_eval, f1.eval(DR_x) + f2.eval(DR_x)]; %record values of objective function - DR_y = DR_y + paramsolver.lambda*(f1.prox(2*DR_x-DR_y, paramsolver.gamma)-DR_x); - if cnt > 1 - relat_change_coefs = norm(DR_x-DR_x_old) / norm(DR_x_old); - relat_change_obj = norm(obj_eval(end) - obj_eval(end-1)) / norm(obj_eval(end-1)); - if paramsolver.verbose > 1 - fprintf(' relative change in coefficients: %e \n', relat_change_coefs); - fprintf(' relative change in objective fun: %e \n', relat_change_obj); - fprintf('\n'); - end - end - DR_x_old = DR_x; - cnt = cnt + 1; - - end - - DR_x = f2.prox(DR_y); %final projection into the constraints - y_dequant = A(DR_x); %dequantized signal - - disp(['Finished after ' num2str(cnt) ' iterations.']) - - - %compare UNLOCBOX with manual - figure - plot([y_dequant A(sol)]) - norm(y_dequant - A(sol)) - title('UNLOCBOX vs. manual solution') + y_dequant = A(sol); - %behaviour of objective through iterations - figure - plot(obj_eval) - title('Objective function value (after projection into constraints)') +% pause +% +% %%%%%%%% manual version %%%%%%%%%%%%% +% %starting point +% DR_y = A(y_quant); +% DR_x_old = DR_y; +% +% relat_change_coefs = 1; +% relat_change_obj = 1; +% cnt = 1; %iteration counter +% obj_eval = []; +% +% % while relat_change_coefs > 0.00001 +% while relat_change_obj > paramsolver.tol +% % DR: gamma = 1 +% DR_x = f2.prox(DR_y,[]); +% obj_eval = [obj_eval, f1.eval(DR_x) + f2.eval(DR_x)]; %record values of objective function +% DR_y = DR_y + paramsolver.lambda*(f1.prox(2*DR_x-DR_y, paramsolver.gamma)-DR_x); +% if cnt > 1 +% relat_change_coefs = norm(DR_x-DR_x_old) / norm(DR_x_old); +% relat_change_obj = norm(obj_eval(end) - obj_eval(end-1)) / norm(obj_eval(end-1)); +% if paramsolver.verbose > 1 +% fprintf(' relative change in coefficients: %e \n', relat_change_coefs); +% fprintf(' relative change in objective fun: %e \n', relat_change_obj); +% fprintf('\n'); +% end +% end +% DR_x_old = DR_x; +% cnt = cnt + 1; +% +% end +% +% DR_x = f2.prox(DR_y); %final projection into the constraints +% y_dequant = A(DR_x); %dequantized signal +% +% disp(['Finished after ' num2str(cnt) ' iterations.']) +% +% +% %compare UNLOCBOX with manual +% figure +% plot([y_dequant A(sol)]) +% norm(y_dequant - A(sol)) +% title('UNLOCBOX vs. manual solution') +% +% %behaviour of objective through iterations +% figure +% plot(obj_eval) +% title('Objective function value (after projection into constraints in each iteration)') end @@ -242,7 +283,7 @@ axis tight legend([h_quant_noise h_dequant_error], 'Quantizat. error', 'Error of reconstr.') -%coefficients of recontructed signal +%coefficients of reconstructed signal figure(fig_coef) hold on h_coefs_dequant = bar(sol,'FaceColor',green);