thfig26.m

％ Atomizer Main Directory, Version .802 里面信号含有分解去噪合成过程的代码 ％---------------------------------------

% thfig26: BP Thesis Figure 2.6 -- Conter examples for MP and OMP
% Signal Length: 1024
% Observations:
% 	(a) DeVore&Temlyakov's greedy death for MP: compression number plot
%	
% 	Shows that rate of convergence need be no faster
% 	than k^{-1/2} in general, even on objects
% 	whose optimal representations contain only m=10 terms.
%
% 	The dictionary is the Dirac dictionary plus a single term, g.
% 	The smallest nonzero singular value of the dictionary
% 	is 1.0. The Signal f is the equal weighted sum of the first
%	ten diracs.
%
% 	The idea is that g is a bit closer to f than either of
% 	the ten vectors in the ideal decomposition.  Once a
% 	mistake is made, the method is completely off the mark;
% 	it keeps going until it takes in the whole dictionary.
%
% 	(b) Shaobing Chen's example for OMP: compression number plot
%
%	Shows that rate of convergence is k^{-1/2}, even on objects
% 	whose optimal representations contain only m=10 terms.
%
%	The dictionary is complete. The first 10 atoms of the dictionary
%	come from the Dirac dictionary. The signal is a simple
%	equi-weighted linear combination of the  first 10 Dirac atoms. The
%	rest n-10 atoms of the dictionary are linear combination of the
%	correponding Dirac and the signal.
%
%	OMF totally fails on Chen's example. It choose all atoms except
%	the first 10 before ever choosing one of the first 10.
%

help thfig26.m

n=1024; m=10;
f = zeros(n,1);
f(1:m) = ones(m,1); f = f / norm(f);
g = zeros(n,1);
g(m:n) = (1 ./ (1:(n-m+1)));
xi        = norm(g);
g  = sqrt(m^2 - m - .05) .* g ./ xi;
g(1:m) = ones(m,1);
g = g ./ norm(g);
residual = f;
history = zeros(n,4);
for i=1:1024,
        history(i,1) = norm(residual);
        [maxcomp,maxind] = max(abs(residual));
        gcoef = g' * residual; gsize = abs(gcoef);
        %plot([ gcoef ; abs(residual(1:25))]); drawnow;
        history(i,1) = norm(residual);
        history(i,2) = gsize;
        history(i,3) = maxcomp;
        history(i,4) = maxind(1);
%       fprintf('step(%i), norm(%g),',i,history(i,1))
        if maxcomp > gsize,
                ind   = maxind(1);
                residual(ind) = 0;
%               fprintf('remove(%i)\n',ind)
        else
                residual = residual - gcoef .* g;
%               fprintf('remove(g)\n')
        end
                
end

figure(1);clf
subplot(2,1,1);
loglog(history(1:1024,1), '--')
ideal = [(9:-1:1) / m 1e-30];
axis([1 1e4 1e-4 1])
hold on;
loglog(ideal);
title('(a) MP on DeVore and Temlyakov''s example')
xlabel('m, Number of Terms in Reconstruction')
ylabel('Reconstruction Error')
axes('Position',[.75 .65 .9 .9],'Visible','off');
h1 = text(0, 0, 'Ideal: Solid');
set(h1, 'FontSize', 6);
h2 = text(0, .05, 'Greedy: Dashed');
set(h2, 'FontSize', 6);
hold off;



%%Chen's example for OMP
n = 1024; m = 10;
theta = [ones(1,m) * 1 / m^.5 zeros(1,n-m)];

a = 2/m;
b = 1 - a;
a = a .^ .5;
b = b .^ .5;

resOMP = zeros(n,1);
for i = 1:n,
	resOMP(i) = ((1 - a^2) / (1 - a^2 + i * a^2)) ^ .5;
end

subplot(2,1,2);
loglog(resOMP, '--');
title('(b) OMP on Chen''s example')
xlabel('m, Number of Terms in Reconstruction')
ylabel('Reconstruction Error')
ideal = [(9:-1:1) / m 1e-30];
axis([1 1e4 1e-4 1])
hold on;
loglog(ideal);
axes('Position',[.75 .15 .9 .9],'Visible','off');
h1 = text(0, 0, 'Ideal: Solid');
set(h1, 'FontSize', 6);
h2 = text(0, .05, 'Greedy: Dashed');
set(h2, 'FontSize', 6);
hold off;