examp.m

这是在网上下的一个东东

%
% Reset to a known state.
%
clear
rand('state',0);
randn('state',0);
%
% setup our problem and target model.
%
[G,b,x]=shaw(20);
mtarget=ones(20,1);
mtarget(1:10)=2*ones(10,1);
%
% Plot the target model.
%
figure(1);
bookfonts;
plotconst(mtarget,-pi/2,pi/2);
axis([-2 2 0 3]);
xlabel('\theta');
ylabel('Intensity');
bookfonts;
%print -deps mtarget.eps
%
% Get the data and plot it.
%
d=G*mtarget;
dn=d+0.001*randn(20,1)
figure(2);
bookfonts;
plotconst(dn,-pi/2,pi/2);
axis([-2 2 0 7]);
xlabel('\theta');
ylabel('Intensity');
%print -deps mdata.eps
%
% Do zeroth order Tikhonov regularization
%
[U,S,V]=svd(G);
figure(3);
bookfonts;
reg_corner=l_curve(U,diag(S),dn,'Tikh');
mtik0=tikhonov(U,diag(S),V,dn,reg_corner);
figure(4);
bookfonts;
plotconst(mtik0,-pi/2,pi/2);
axis([-2 2 0 3]);
xlabel('\theta');
ylabel('Intensity');
%print -deps ppmtik0.eps
%
% Do second order regularization.
%
L=full(get_l(20,2));
[U,sm,X,V]=cgsvd(G,L);
figure(5);
bookfonts;
reg_corner2=l_curve(U,sm,dn,'Tikh');
mtik2=tikhonov(U,sm,X,dn,reg_corner2);
figure(6);
bookfonts;
plotconst(mtik2,-pi/2,pi/2);
axis([-2 2 0 3]);
xlabel('\theta');
ylabel('Intensity');
%print -deps ppmtik2.eps
%
% Compute a TV regularized solution with IRLS.  
%
alpha=1.0;
L1=full(get_l(20,1));
A=[G; alpha*L1];
b=[d; zeros(19,1)];
figure(7);
bookfonts;
mtv=irls(A,b,1.0e-3,1.0e-5,1,100);
plotconst(mtv,-pi/2,pi/2);
axis([-2 2 0 3]);
xlabel('\theta');
ylabel('Intensity');
%print -deps mtv.eps
disp('mtv quality');
norm(G*mtv-dn)
norm(mtv-mtarget)/norm(mtarget)
%
% Compute the pptsvd solution.
%
mpptsvd=pptsvds(G,L1,dn,2);
figure(8);
bookfonts;
plotconst(mpptsvd,-pi/2,pi/2);
axis([-2 2 0 3]);
xlabel('\theta');
ylabel('Intensity');
%print -deps mpptsvd.eps
disp('mpptsvd quality');
norm(G*mpptsvd-dn)
norm(mpptsvd-mtarget)/norm(mtarget)