examp.m

这是在网上下的一个东东

%
% This script generates the Shaw example from chapter 4 of the book.
%
%
% First, clear and reset the RNG states.
%
clear
randn('state',0);
rand('state',0);
%
% Build the n=20 Shaw problem.
%
[G,m,d]=shaw(20);
%
% Compute the SVD.
%
[U,S,V]=svd(G);
%
% Plot the singular values of G.
%
clf;
bookfonts
semilogy(diag(S),'ko');
xlabel('i');
ylabel('s_i');
%print -deps c4fspec20fig.eps
%
% Plot column 18 of V.
%
figure(2);
clf;
bookfonts
plotconst(V(:,18),-pi/2,pi/2);
xlabel('\theta (radians)')
ylabel('Intensity')
%print -deps shawv18.eps
%
% plot column 1 of V.
%
figure(3);
clf;
bookfonts
plotconst(V(:,1),-pi/2,pi/2);
xlabel('\theta (radians)')
ylabel('Intensity')
%print -deps shawv1.eps
%
% plot the spike model.
%
spike=zeros(20,1);
spike(10)=1;
figure(4);
clf;
bookfonts
plotconst(spike,-pi/2,pi/2);
ylabel('Intensity')
xlabel('\theta (radians)')
axis([-2 2 -0.5 1.5]);
%print -deps shawspikemod.eps
%
% Plot the shawspike data.
%
spikedata=G*spike;
%
% Chpater 5 script uses dspike, not spikedata.
%
dspike=spikedata;
figure(5);
clf;
bookfonts
plotconst(spikedata,-pi/2,pi/2);
ylabel('Intensity')
xlabel('s (radians)')
axis([-2 2 -0.5 1.5]);
%print -deps shawspikedif.eps
%
% Create slightly noisy data and see what happens.
%
dn=spikedata+1.0e-6*randn(size(spikedata));
%
% Because the chapter 5 example script uses dspiken rather than dn...
%
%solution for no noise
spikemodnn=G\spikedata;
dspiken=dn;
%
% Find the pseudoinverse solution, p=18.
%
p=18;
Up=U(:,1:p);
Vp=V(:,1:p);
Sp=S(1:p,1:p);
Gg=Vp*inv(Sp)*Up';
spikemod18=Gg*dn;
%
%
figure(6);
clf;
bookfonts
plotconst(spikemodnn,-pi/2,pi/2);
ylabel('Intensity')
xlabel('\theta (radians)')
%print -deps shawspikeinvnn.eps

figure(7);
clf;
bookfonts
plotconst(spikemod18,-pi/2,pi/2);
ylabel('Intensity')
xlabel('\theta (radians)')
%print -deps shawspikeinv18.eps

%
% Now, p=10.
%
p=10;
Up=U(:,1:p);
Vp=V(:,1:p);
Sp=S(1:p,1:p);
Gg=Vp*inv(Sp)*Up';
spikemod10=Gg*spikedata;
%
%
figure(8);
clf;
bookfonts
plotconst(spikemod10,-pi/2,pi/2);
ylabel('Intensity')
xlabel('\theta (radians)')
%print -deps shawspikeinv10nn.eps
%
% Now, the noisy version.
%
spikemod10n=Gg*dn;
figure(9);
clf;
bookfonts
plotconst(spikemod10n,-pi/2,pi/2);
ylabel('Intensity')
xlabel('\theta (radians)')
%print -deps shawspikeinv10.eps
%
% Plot the spectrum for n=100 problem.
%
[G100,m100,d100]=shaw(100);
[U100,S100,V100]=svd(G100);
figure(10);
clf;
bookfonts
semilogy(diag(S100),'ko');
xlabel('i')
ylabel('s_i')
%print -deps shawspec100.eps
%
% 
%
spike100=zeros(100,1);
spike100(46:50)=ones(5,1);
spikedata100=G100*spike100;
spikedata100n=spikedata100+1.0e-6*randn(100,1);
p=10;
Up=U100(:,1:p);
Vp=V100(:,1:p);
Sp=S100(1:p,1:p);
Gg=Vp*inv(Sp)*Up';
spikeinv100n=Gg*spikedata100n;
figure(11);
clf;
bookfonts
plotconst(spikeinv100n,-pi/2,pi/2);
ylabel('Intensity')
xlabel('\theta (radians)')
%print -deps spikeinv100.eps
%
%
%
p=18;
Up=U100(:,1:p);
Vp=V100(:,1:p);
Sp=S100(1:p,1:p);
Gg=Vp*inv(Sp)*Up';
spikeinv100n18=Gg*spikedata100n;
figure(12);
clf;
bookfonts
plotconst(spikeinv100n18,-pi/2,pi/2);
ylabel('Intensity')
xlabel('\theta (radians)')
%print -deps spikeinv100n18.eps
%
% Now, the 6 element discretization.
%
[G6,m6,d6]=shaw(6);
[U6,S6,V6]=svd(G6);
figure(13);
clf;
bookfonts
semilogy(diag(S6),'ko');
xlabel('i')
ylabel('s_i')
%print -deps shawspec6.eps
%
% Save this for use in later examples.
%
%save shawexamp.mat