function svd_demo (n,p)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SVD_DEMO
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for MTH 451-551
%% this demo shows an application of SVD (Singular Value Decomposition) 
%% to image processing,
%% references: [J. Daniel, B. Noble, "Applied Linear Algebra"]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% M. Peszynska,  2004/11
%% Copyright Department of Mathematics, Oregon State University
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% INSTRUCTIONS: 
%% calling sequence: 
%%        svd_demo(n,p)
%% n x n is the size of the image 
%% p is the thickness of the band
%% default values are n=20,p=1
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
if (nargin ~= 2) 
   n = 20;          %% default size of the image
   p = 1;           %% default thickness of the band
end;

%%%%%%%%% first step: create a discretized "image":
%%             here we discretize the image of a letter 'x'
a = diag(ones(n,1),0);
for i = 1:p s = ones(n-i,1); a = a + diag(s,i) + diag(s,-i);end;
for i = 1:n/2 for j = 1:n/2 a(n-j+1,i) = a(j,i); end;end;
for i = 1:n/2+1 for j = 1:n/2+1 a(i,n-j+1) = a(i,j); end;end;
%% display the matrix
spy(a); pause;
pcolor(a); pause;

%%%%%%%%% find SVD of a
[u,s,v] = svd(a);
diag(s),                %% show the singular values of a
m = rank(a);              %% what is the rank of a ?

%%%%%%%%% find successive low rank approximations of a in matrix <appr>
appr = zeros(n,n);
for i=1:m 
    appr=appr+u(:,i)*s(i,i)*v(:,i)';
                          %% report how good is the current approximation
    'approximation of order ',i,' error ',norm(a-appr,2)/norm(a),
    pcolor(appr);pause,   %% and display it

end;