omp_matrix.m

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

function coef = OMP_Matrix(x, F, natom, frac)
% OMP_Matrix -- Orthogonal Matching Pursuit with direct matrix-on-vector 
%		operations for any dictionary
%  Usage
%	coef = OMP_Matrix(x, F[, natom, frac])
%  Input
%	x	1-d signal; n-by-1 column vector
%	F	dictionary matrix; each columns corresponds to each atom
%	natom   max # of atoms desired, default = n
%	frac    min fraction total signal energy to enter, default=1e-2
%  Outputs
%	coef	coef of the MP representation; column vector
%  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, OMP, MP_Matrix
%

x = x(:);
[n, N] = size(F);
if nargin < 4,
	natom = N;
	frac = 1e-2;
end
res = [];
index = [];
GramInvNew = zeros(natom);
k = 0;
residule = x;
nrm = norm(x);
amp = nrm;
fprintf('\nMPBF_Matrix:\n');
while (amp > frac*nrm) & (k < natom),
	P = abs(F' * residule);
	[amp i] = max(P);
	i = i(1);
	newterm = F(:, i);
	index = [index i];
	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 = F(:, index);
	residule = x - X * GramInv * X' * x;
	resnorm = norm(residule);
	res = [res resnorm];
	disp(sprintf('Step%4g:	select = %g	norm(residule) = %g', k, i, resnorm));
	if resnorm < 1e-10, break; end

end

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

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