bpfig25.m

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

% bpfig25: BP Figure 2.5 -- Counter-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 remaining n-10 atoms of the dictionary are a linear
%	combination of the correponding Dirac and the signal.
%
%	OMP totally fails on Chen's example. It choose all atoms except
%	the first 10 before ever choosing one of the first 10.
%----------------------------------------------------------------------

%----------------------------------------------------------------------
%        1998: (S. Chen) Original script for BP paper.
% 09 Apr 2001: (M. Saunders) Minor editing.
%----------------------------------------------------------------------

help bpfig25
bpengine;

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      = g * (sqrt(m^2 - m - .05) / 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

%--------------------
% Plots
%--------------------
fprintf('\nConstructing figure(%1g) ...\n', FIGURE)
figure(FIGURE);   clf reset;

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) / 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;