mof.m

来自「％ Atomizer Main Directory, Version .802 」· M 代码 · 共 72 行

function c = MOF(x, NameOfDict, par1, par2, par3, CGAccuracy)
% MOF.m -- Method of Frames using a CG solver:
%		Min || c ||_2  subject to \Phi * c = x
%	i.e.
%		c = \Phi^T * (\Phi * \Phi^T)^ (-1) * x
%
%  Usage
%	c = MOF(x, NameOfDict, par1, par2, par3[, CGAccuracy])
%  Inputs
%	x		1-d signal; n-by-1 column vector
%	NameOfDict	string; name of the dictionary
%	par1,par2,par3	parameters of the dictionary
%	CGAccuracy	accuracy of the CG solver, default = 1e-5
%
%	Use 'help dictionary' for dictionary objects: NameOfDict,par1,par2,par3
%  Description
%	1. If the dictionary is a tight frame, then CG converges within
%	   one step and achieves accuracy 1e-16.
%	2. Paramter CGAccuracy controls the accuracy of the CG solver, and
%	   hence the accuracy of MOF. One may want to reset its value 
%	   according to his purpose.
%	3. Paramter MaxCGIter controls the maximum number of CG iterations.
%	   The default m(the cadinality of the dictionary) is often more 
%	   than enough. But if one may want to increase this parameter 
%	   when he gets message saying "Too many CG iterations".
%  Outputs
%	c		coef of the MOF representation; column vector
%  See Also
%	MOF_Matrix
%

	%% Dimensions of the problem
	x = x(:);
	n = length(x);
	[m L] = SizeOfDict(n, NameOfDict, par1, par2, par3);

	%% CG parameters
	if nargin < 6,
		CGAccuracy = 1e-5;
	end
	MaxCGIter = m;

	%% Conjugate Gradient Method for (Phi*Phi') y = x
	%fprintf('\nMOF:\n');
	thresh = CGAccuracy * norm(x); thresh = thresh ^ 2;
	y = zeros(n,1);
	r = x;
	rho_1 = norm(r) ^ 2;
	for k = 1:MaxCGIter,
		%fprintf('	CG Step: %3g:	%g\n', k, rho_1);
		if k == 1,
			p = r;
		else
			beta = rho_1 / rho_0;
			p = r + beta * p;
		end
		w = FastS(FastA(p, NameOfDict, par1, par2, par3), n, ...
			NameOfDict, par1, par2, par3);
		alpha = rho_1 / (p' * w);
		r = r - alpha * w;
		rho_0 = rho_1;
		rho_1 = norm(r) ^2;
		y = y + alpha * p;
		if rho_1 < thresh,
			break;
		end
	end
	if k == MaxCGIter,
		disp('Too many CG iterations!')
	end
	c = FastA(y, NameOfDict, par1, par2, par3);