[c3bd51]: frames / frsyniter.m  Maximize  Restore  History

Download this file

144 lines (118 with data), 4.2 kB

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
function [f,relres,iter]=frsyniter(F,c,varargin)
%FRSYNITER Iterative synthesis
% Usage: f=frsyniter(F,c);
%
% Input parameters:
% F : Frame
% c : Array of coefficients.
% Ls : length of signal.
% Output parameters:
% f : Signal.
% relres : Vector of residuals.
% iter : Number of iterations done.
%
% `f=frsyniter(F,c)` iteratively inverts the analysis operator of *F*, so
% `frsyniter` always performs the inverse operation of |frana|, even
% when a perfect reconstruction is not possible by using |frsyn|.
%
% `[f,relres,iter]=frsyniter(...)` additionally returns the relative
% residuals in a vector *relres* and the number of iteration steps *iter*.
%
% **Note:** If it is possible to explicitly calculate the canonical dual
% frame then this is usually a much faster method than invoking
% `frsyniter`.
%
% `frsyniter` takes the following parameters at the end of the line of
% input arguments:
%
% 'tol',t Stop if relative residual error is less than the
% specified tolerance. Default is 1e-9
%
% 'maxit',n Do at most n iterations.
%
% 'cg' Solve the problem using the Conjugate Gradient
% algorithm. This is the default.
%
% 'pcg' Solve the problem using the Preconditioned Conjugate Gradient
% algorithm.
%
% 'print' Display the progress.
%
% 'quiet' Don't print anything, this is the default.
%
% Examples
% --------
%
% The following example shows how to rectruct a signal without ever
% using the dual frame:::
%
% F=frame('dgtreal','gauss',10,20);
% c=frana(F,bat);
% [r,relres]=frsyniter(F,c,'tol',1e-14);
% norm(bat-r)/norm(bat)
% semilogy(relres);
% title('Conversion rate of the CG algorithm');
% xlabel('No. of iterations');
% ylabel('Relative residual');
%
% See also: frame, frana, frsyn, franaiter
% AUTHORS: Nicki Holighaus & Peter L. S��ndergaard
if nargin<2
error('%s: Too few input parameters.',upper(mfilename));
end;
definput.keyvals.Ls=[];
definput.keyvals.tol=1e-9;
definput.keyvals.maxit=100;
definput.flags.alg={'cg','pcg'};
definput.keyvals.printstep=10;
definput.flags.print={'quiet','print'};
[flags,kv,Ls]=ltfatarghelper({'Ls'},definput,varargin);
% Determine L from the first vector, it must match for all of them.
L=framelengthcoef(F,size(c,1));
A=@(x) frsyn(F,frana(F,x));
% It is possible to specify the initial guess, but this is not
% currently done
if flags.do_pcg
d=framediag(F,L);
M=spdiags(d,0,L,L);
[f,flag,~,iter,relres]=pcg(A,frsyn(F,c),kv.tol,kv.maxit,M);
else
[f,flag,~,iter,relres]=pcg(A,frsyn(F,c),kv.tol,kv.maxit);
end;
if nargout>1
relres=relres/norm(c(:));
end;
% Cut or extend f to the correct length, if desired.
if ~isempty(Ls)
f=postpad(f,Ls);
else
Ls=L;
end;
if 0
% This code has been disabled, as the PCG algorithm is so much faster.
if flags.do_unlocbox
% Get the upper frame bound (Or an estimation bigger than the bound)
[~,B]=framebounds(F,L,'a');
% Set the parameter for the fast projection on a B2 ball
param.At=@(x) frsyn(F,x); % adjoint operator
param.A=@(x) frana(F,x); % direct operator
param.y=c; % coefficient
param.tight=0; % It's not a tight frame
param.max_iter=kv.maxit;
param.tol=kv.tol;
param.nu=B;
% Display parameter 0 nothing, 1 summary at convergence, 2 all
% steps
if flags.do_print
param.verbose=1;
else
param.verbose=0;
end;
% Make the projection. Requires UNLocBOX
[f, ~] = fast_proj_B2(zeros(L,1), 0, param);
% compute the residue
res = param.A(f) - param.y; norm_res = norm(res(:), 2);
relres=norm_res/norm(c(:), 2);
iter=0; % The code of the fast_proj_B2 is not yet compatible with this
end;
end;

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:

JavaScript is required for this form.





No, thanks