%%
%% Code for demos presented in lecture 2
%% Slightly improved after the lecture
%%
%% P. Perona and M. Weber - April 1 1999 (not a joke!)
%% EE/CNS 148 - Spring 1999
%% (c) California Institute of Technology
%% 

%%
%%create image
%%

%% Standard deviation of noise:
NOISE_STD = 1;  %%% <- Change this value in order to test the robustness of the method

%% A few formatting parameters and initialization of ima array
LN_WDTH = 2; 	%% Thickness of line in plots
IMA_LEN = 128; %% Size of the 1D image
ima = zeros(1,IMA_LEN);

%% Generate signal
s = [-1.1*ones(1,10) 2.1*ones(1,6)]; %% Some signal, you may change this
SIG_LEN = length(s);
SIG_POS = IMA_LEN/2;
ima(SIG_POS+[1:SIG_LEN]) = s; %% Insert the signal into the image array
figure(1); subplot(3,1,1); plot(ima,'LineWidth',LN_WDTH);
title('Image with signal and no noise.');

%% Add noise to image
ima = ima + NOISE_STD*randn(1,IMA_LEN); %% This is white Gaussian noise
figure(1); subplot(3,1,2); plot(ima,'LineWidth',LN_WDTH);
title(['Image with white Gaussian noise of std dev \sigma=' num2str(NOISE_STD)]);


%%
%% Match-filter the image and calculate an appropriate threshold
%%
k = fliplr(s);           %% Set the kernel to be the matched filter
                         %% It needs to be flipped because of how conv2() will use it
sig_norm = sqrt(sum(k.^2));
k = k / sig_norm; %% Normalize the kernek k
filter_out = conv2(ima,k,'same');
figure(1); subplot(3,1,3); plot(filter_out,'LineWidth',LN_WDTH);
THRESH = (sig_norm - 0)/2;
hold on; plot([1 IMA_LEN],[THRESH THRESH],'g','LineWidth',LN_WDTH); hold off;
title('Output of matched filter; the thresh is set at half-point between the E d and E d_n');



%%%%%%%%%%%%%
%%%%%%%%%%%%%
%%
%% Variable shape signal experiment - principal component analysis
%%
%%%%%%%%%%%%%
%%%%%%%%%%%%%
ima = zeros(1,IMA_LEN);
MAX_SIG_LEN = 32;
MAX_JITTER = 8;
N = 100;
S = zeros(N,MAX_SIG_LEN);

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

figure(5); title('16 sample signals');
for i=1:16,
   subplot(4,4,i);
   plot(S(i,:));
end;


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

figure(6); sval = diag(SS); semilogy(sval(1:10),'.','MarkerSize',15);
title('The first singular values of S S^T on a log scale');
   
figure(7); title('The first principal components (eigenvectors)');
for i=1:16,
   subplot(4,4,i);
   plot(U(:,i));
   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);
end;
figure(8); plot(c1,c2,'.','MarkerSize',15);
xlabel('1st component'); ylabel('2nd component');
title(['The ' num2str(N) ' sample signals projected onto the reduced base']);
axis('equal'); axis([-10 10 -10 10]);