Homec4science
Diffusion iCAPs public version 2 (master) Browse

iCAPs_v2/packages/gspbox/prox/gsp_wavelet_dn.m
e5f5005bb2efmaster

Last Change
Raw File

gsp_wavelet_dn.m
View Options

function [sol, info] = gsp_wavelet_dn(G, w, x, lambda, param)
%GSP_WAVELET_DN Wavelet denoising
% Usage: sol = gsp_wavelet_dn(G, w, x, lambda, param)
% sol = gsp_wavelet_dn(G, w, x, lambda)
% [sol, info] = gsp_wavelet_dn(...)
%
% Input parameters:
% G : Graph structure
% w : Wavelet filterbank
% x : Signal to be denoised
% lambda: Regularization parameter
% param : Structure of optional parameters.
% Output parameters
% sol : Solution.
% info : Structure summarizing informations at convergence
%
% This function will denoise a signal by solving the following convex
% promblem:
%
% sol = argmin_{z} || W^* z - x' ||_2^2 + lambda * ||z||_1
%
% Where W is the frame associated to the filterbank w, x' a part of
% the signal to be denoised and z the wavelet coefficient.
%
% x' consists of the high frequency part of x. It is obtained by
% setting down to zero the low pass filter of the filerbank w.
%
% This function require the UNLocBoX to be executed. You can download it
% at http://unlcobox.sourceforge.net
%
% param is a Matlab structure containing the following fields:
%
% param.tol : is stop criterion for the loop. The algorithm stops if
%
% ( n(t) - n(t-1) ) / n(t) < tol,
%
% where n(t) = f(x)+ 0.5 X-Z_2^2 is the objective function at iteration t*
% by default, tol=10e-4.
%
% param.maxit : max. nb. of iterations (default: 200).
%
% param.tight : 1 W^ are both tight frame or 0 if not
% (default = 0)
%
% param.verbose : 0 no log, 1 a summary at convergence, 2 print main
% steps (default: 1)
%
% param.method : Solver to be used ('FISTA', 'ISTA', 'DG') By default
% it is 'FISTA'. ('DG' is Douglas Rachford)
%
% param.gamma : stepsize for the 'DG' algorithm
%
% info is a Matlab structure containing the following fields:
%
% info.algo : Algorithm used
%
% info.iter : Number of iteration
%
% info.time : Time of exectution of the function in sec.
%
% info.final_eval : Final evaluation of the function
%
% info.crit : Stopping critterion used
%
% Demo: gsp_demo_wavelet_dn
%
% See also: gsp_prox_tv gsp_solve_l1 gsp_solve_l0
%
%
% Url: http://lts2research.epfl.ch/gsp/doc/prox/gsp_wavelet_dn.php
% Copyright (C) 2013-2016 Nathanael Perraudin, Johan Paratte, David I Shuman.
% This file is part of GSPbox version 0.7.0
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
% If you use this toolbox please kindly cite
% N. Perraudin, J. Paratte, D. Shuman, V. Kalofolias, P. Vandergheynst,
% and D. K. Hammond. GSPBOX: A toolbox for signal processing on graphs.
% ArXiv e-prints, Aug. 2014.
% http://arxiv.org/abs/1408.5781
% Author: Nathanael Perraudin
% Date: 25 March 2014
% Testing:
if nargin<5
param = struct;
end
if ~isfield(param, 'verbose'), param.verbose = 2; end
if ~isfield(param, 'tight'), param.tight = 0; end
if ~isfield(param, 'tol'), param.tol = 1e-3; end
if ~isfield(param, 'maxit'), param.maxit = 200; end
if ~isfield(param, 'method'), param.method = 'FISTA'; end
if ~isfield(param, 'gamma'), param.gamma = 1; end
if ~isfield(param, 'normalize'), param.normalize = 0; end
if param.normalize
ntig = gsp_norm_tig(G,w(2:end));
ntig = gsp_mat2vec(ntig);
else
ntig = 1;
end
param_solver.maxit = param.maxit;
param_solver.verbose = param.verbose;
param_solver.tol = param.tol;
% % First method to handle low frequencies
% N = G.N;
% % Remove the low frequency part
% % 1) Analysis
% y_noise_c = gsp_filter_analysis(G,w,x);
% % % 2) Save low frequencies
% % y_noise_c_bf = zeros(size(y_noise_c));
% % y_noise_c_bf(1:N) = y_noise_c(1:N);
% % 3) Set low frequencies to 0
% y_noise_c(1:N) = 0;
% % 4) Inverse to find the signal without low frequencies
% y_noise_2 = gsp_filter_inverse(G,w,y_noise_c);
% y_noise_bf = x - y_noise_2;
% Remove the low frequency part
lf = gsp_design_expwin(G,0.3);
% 1) Analysis
y_noise_bf = gsp_filter_analysis(G,lf,x);
% 2) Save low frequencies
y_noise_2 = x-y_noise_bf;
% Sparsity of the coefficients -- prox l1
param_l1.verbose = param.verbose -1;
fc1.prox = @(x,T) prox_l1(x,T*lambda,param_l1);
fc1.eval = @(x) lambda*norm(x(:),1);
% Sparsity on the coefficients side -- prox l2
[~, B] = gsp_filterbank_bounds(G,w);
param_l2f.verbose = param.verbose - 1;
param_l2f.y = y_noise_2;
param_l2f.nu = 1/(2*B^2);
param_l2f.tight = param.tight;
param_l2f.ntig = ntig;
Wh = w(2:end); % drop the low pass filter
fl2.prox = @(x,T) gsp_prox_l2_filterbank(x, T, G, Wh, param_l2f);
fl2.eval = @(x) gsp_norm_l2_filterbank(G,Wh,x,param_l2f);
Wop = @(x) gsp_filter_analysis(G, Wh, x)./ntig;
Wtop = @(x) gsp_filter_synthesis(G, Wh, x./ntig);
fl2.grad = @(x) 2*Wop(Wtop(x)-y_noise_2);
fl2.beta = (2*B^2);
% starting point
x0 = Wop(x);
switch param.method
case 'FISTA'
param_solver.gamma = min(1/(2*B^2),0.5)*(min(ntig)^2);
[denoise_y_wav2_c,info] = forward_backward(x0,fc1,fl2,param_solver);
case 'ISTA'
param_solver.gamma = min(1/(2*B^2),0.5)*(min(ntig)^2);
param_solver.method = 'ISTA';
[denoise_y_wav2_c,info] = forward_backward(x0,fc1,fl2,param_solver);
case 'DG'
param_solver.gamma = param.gamma;
[denoise_y_wav2_c,info] = douglas_rachford(x0,fc1,fl2,param_solver);
otherwise
error('GSP_WAVELET_DN: Unknown method! ');
end
sol =Wtop(denoise_y_wav2_c)+y_noise_bf;
end
c4science · Help