omp.m

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

function coef = OMP(x, NameOfDict, par1, par2, par3, natom, frac)
% OMP -- Orthogonal Matching Pursuit
%  Usage
%	coef = OMP(x, NameOfDict, par1, par2, par3[, natom, frac])
% Input
%	x		1-d signal; n-by-1 column vector
%	NameOfDict	string; name of the dictionary
%	par1,par2,par3	parameters of the dictionary
%	natom    	max # of atoms desired, default = 100
%	frac     	min fraction total signal energy to enter, default=1e-2
%
%	Use 'help dictionary' for dictionary objects: NameOfDict,par1,par2,par3
%  Outputs
%	coef		coef of the MP representation
%  Description
%	1. OMP is like MP, except each step is followed by a backfitting
%	   (least square) procedure
%	2. natom controls the maximum number of atoms MP can select
%	3. Selected atoms must have coefficients greater than a certain
%	   fraction of total signal energy, defined by the parameter
%	   frac, in order to be enter.
%  See Also
%	MP, MP_Matrix, OMP_Matrix
%

x = x(:);
n = length(x);

if nargin < 7,
	frac = 1e-2;
end
if nargin < 6,
	natom = 100;
end


[m L] = SizeOfDict(n, NameOfDict, par1, par2, par3);
zerosm = zeros(m,1);
nrm = norm(x);
amp = nrm;

residule = x;
coef = [];
index = [];
k = 0;
GramInvNew = zeros(natom);
X = [];
fprintf('\nOMP:\n');
while (amp > frac*nrm) & (k < natom),
	c = FastA(residule, NameOfDict, par1, par2, par3);
	[amp i] = max(abs(c));
	i = i(1);
	index = [index i];
	temp = zerosm; temp(i) = 1;
	newterm = FastS(temp, n, NameOfDict, par1, par2, par3);
	if k == 0,
		%Initializing everything
		GramInvNew(1,1) = 1;
		k = 1;
		GramInv = 1;
	else
		%Update the Inverse of the Gram Matrix
		Q = GramInv * X' * newterm;
		m = 1 / (newterm' * newterm - ...
			 newterm' * X * GramInv * X' * newterm);
		GramInvNew(1:k, 1:k) = GramInv + Q * Q' * m;
		GramInvNew(1:k, k+1) =  - Q * m;
		GramInvNew(k+1, 1:k) =  - Q' * m;
		GramInvNew(k+1, k+1) = m;
		k = k + 1;
		GramInv = GramInvNew(1:k , 1:k);
	end	

	X = [X newterm];
	residule = x - X * GramInv * X' * x;
	resnorm = norm(residule);
	disp(sprintf('Step%4g:	select = %5g     norm(residule) = %10.2e', k, i, resnorm));
end

coef = GramInv * X' * x;
xrec = X * coef;

coefbest = zeros(m, 1);
coefbest(index) = coef;
coef = coefbest;