mofdn.m

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

function [xrec, c] = MOFDN(y, NameOfDict, par1, par2, par3, CGAccuracy, lambda)
% MOFDeNoise -- MOF DeNoising (Ridge Regression) using a CG slover:
%
%		Min || y - \Phi * c ||_2^2 + lambda * || c ||_2^2
%		
%		i.e.
%			c = (lambda * I + \Phi^T * \Phi) ^ (-1) * \Phi * y
%
% Usage:
%	[xrec, c] = MOFDN(y, NameOfDict, par1, par2, par3
%			[, CGAccuracy, lambda])
% Inputs:
%       y               1-d signal to be de-noised; column vector
%	NameOfDict	string; name of the dictionary
%	par1,par2,par3	parameters of the dictionary
%	CGAccuracy	accuracy of the CG solver, default = 1e-5
%	lambda		regularization parameter, default=(2*log(m))^.5
%			where m = the cardinality of the dictionary
%
%	Use 'help dictionary' for dictionary objects: NameOfDict,par1,par2,par3
% Outputs:
%	xrec		cleaned signal; column vector
%	c		coeff of cleaned signal; column vector
% Description
%	1. Assumes noise level == 1
%	2. One may want to use a different regularization paramter lambda.
%	4. Paramter CGAccuracy controls the accuracy of the CG solver and
%	   hence the accuracy of MOF.
%	5. 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".
% See Also
%	MOF
%

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

	%Setting parameters
	if nargin < 8,
		lambda = (2 * log(m)) ^ .5;
	end
	if nargin < 7,
		CGAccuracy = 1e-5;
	end
	MaxCGIter = m;


	%% Conjugate Gradient Method for (lambda*I + \Phi^T*\Phi) c = \Phi^T*y
	c = zeros(m,1);
	r = FastA(y, NameOfDict, par1, par2, par3);
	thresh = CGAccuracy * norm(r); thresh = thresh ^ 2;
	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 = FastA(FastS(p, n, NameOfDict, par1, par2, par3), ...
			NameOfDict, par1, par2, par3) + lambda * p;
		alpha = rho_1 / (p' * w);
		r = r - alpha * w;
		rho_0 = rho_1;
		rho_1 = norm(r) ^2;
		c = c + alpha * p;
		if rho_1  < thresh,
			break;
		end
	end
	if k == MaxCGIter,
		disp('Too many CG iterations!')
	end
	xrec = FastS(c, n, NameOfDict, par1, par2, par3);