%%
%% Code for demos presented in lecture 3
%% (c) Caltech - P. Perona, April 2000
%% Modified April 7, 2003
%%%%%%%%%%%%%
%%%%%%%%%%%%%
%%
%% Variable shape signal experiment - principal component analysis
%%
%%%%%%%%%%%%%
%%%%%%%%%%%%%
clear;

IMA_LEN = 128; %% Size of the 1D image
MAX_SIG_LEN = 32;
MAX_JITTER = [0, 8];
N = 1000;
NOISE_LEVEL = 0.3;

%% Figure out position where figures will be shown
xy_screen = get(0,'ScreenSize'); xy_screen(4) = xy_screen(4);
x0 = xy_screen(1); x1 = (x0-1)+round((xy_screen(3)-xy_screen(1))/2); dx = x1;
y0 = xy_screen(2); y1 = (y0-1)+round((xy_screen(4)-xy_screen(2))/2); dy = y1;
figpos(1,:) = [x0 y1+1 dx dy]; figpos(2,:) = [x1+1 y1+1 dx dy];
figpos(3,:) = [x0 y0 dx dy]; figpos(4,:) = [x1+1 y0 dx dy];

for jj = 1:length(MAX_JITTER),
    fig_offset = (jj-1)*10;
    
    ima = zeros(1,IMA_LEN);
    S = zeros(N,MAX_SIG_LEN);
    
    for EE=1:N,
        %% Generate N realizations of noisy and variable signal signal
        xx = round(MAX_JITTER(jj)*rand(1,3))+[4, (MAX_SIG_LEN/2)-4, MAX_SIG_LEN-8];
        S(EE,xx(1):xx(2))=-1-rand;
        S(EE,xx(2):xx(3)) =1+rand;
        noise = NOISE_LEVEL*randn(1,length(S(EE,:)));
        S(EE,:) = S(EE,:)+noise;
        %S(EE,:) = S(EE,:) - mean(S(EE,:));
    end;
    
    figure(1+fig_offset);
    set(1+fig_offset,'Name','Sample of original Signals');
    for i=1:16,
        subplot(4,4,i);
        plot(S(i,:)); 
        title(num2str(i)); axis off
    end;
    
    [U,SS,V] = svd(S'*S); sval = diag(SS); 
    figure(2+fig_offset); semilogy(sval(1:20),'.','MarkerSize',15);
    title('The first singular values of S S^T on a log scale');
    
    figure(3+fig_offset); 
    set(3+fig_offset,'Name','Principal vectors');
    for i=1:16,
        subplot(4,4,i);
        plot(U(:,i));
        title(num2str(i)); axis off
    end;
    
    %% Projection of signals onto new base
    for EE=1:N,
        c1(EE) = S(EE,:)* U(:,1); c2(EE) = S(EE,:)* U(:,2); %% Project examples of noisy signals
        
        noise = NOISE_LEVEL*randn(1,length(S(EE,:))); %% Generate example of noise with no signal
        c1n(EE) = noise*U(:,1); c2n(EE) = noise*U(:,2); %% Project examples of pure noise
    end;
    
    
    figure(4+fig_offset); clf; plot(c1,c2,'g.','MarkerSize',15);
    hold on; plot(c1n,c2n,'r.','MarkerSize',15); hold off;
    xlabel('1st component'); ylabel('2nd component');
    title(['The ' num2str(N) ' sample signals (green) and noise (red) projected onto the reduced base']);
    axis('equal'); %axis([-10 10 -10 10]); 
    
    
    for ff=1:4, 
        set(ff+fig_offset,'Position',figpos(ff,:),'MenuBar','none'); 
    end;
    pause
    
end;


%%
%%
%%
%%
%%
clear S U V SS c1 c2 c1n c2n
%%
%% Variable shape signal experiment
%%
MAX_SIG_LEN = 32;
MAX_JITTER = 14;
%% MAX_JITTER = 4; %% Try this too
BLIP_LENGTH = 2;
N = 400;
S = zeros(N,MAX_SIG_LEN);

for EE=1:N,
    %% Generate signal
    xx = round(MAX_JITTER*rand(1,2))+[1 MAX_SIG_LEN-MAX_JITTER];
    S(EE,xx(1)+[1:BLIP_LENGTH])=-1-rand;
    S(EE,xx(2)-[1:BLIP_LENGTH])=+1+rand;
    %S(EE,:) = S(EE,:) - mean(S(EE,:));
end;

figure(100+1);
set(100+1, 'Name','Original Signals');
for i=1:16,
    subplot(4,4,i);
    plot(S(i,:)); axis off
    title(num2str(i));
end;

%MS = mean(S);
%S = S - ones(N,1)*MS;


[U,SS,V] = svd(S'*S);

figure(100+2); sval = diag(SS); semilogy(sval,'.','MarkerSize',15);
title('Principal components (singular values)');


figure(100+3);
set(100+3, 'Name','Principal vectors');
for i=1:16,
    subplot(4,4,i);
    plot(U(:,i)); axis off
    title(num2str(i));
end;


%% Projection of signals onto new base
for EE=1:N,
    c1(EE) = S(EE,:) * U(:,1);    c2(EE) = S(EE,:) * U(:,2);
    
    noise = NOISE_LEVEL*randn(1,length(S(EE,:))); %% Generate example of noise with no signal
    c1n(EE) = noise*U(:,1); c2n(EE) = noise*U(:,2); %% Project examples of pure noise
end;
figure(100+4); plot(c1,c2,'g.','MarkerSize',15); hold on;
plot(c1n,c2n,'r.','MarkerSize',15); hold off;
xlabel('1st component'); ylabel('2nd component');
title(['The ' num2str(N) ' sample signals and noise projected onto the reduced base']);
axis equal; axis([-3 3 -3 3]);


for ff=1:4, set(100+ff,'Position',figpos(ff,:),'MenuBar','none'); end;