%%
%% Code for demos presented in lecture 2
%%


%% P. Perona and M. Weber - April 1 1999
%% EE/CNS 148 - Spring 1999
%% Modified March 31 2000
%%
%% (c) California Institute of Technology
%% 

%%
%%create image
%%
%% Standard deviation of noise:
NOISE_STD = [0.1 1 2 4];  %%% <- 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
KSMOOTH = 0;

for i=1:length(NOISE_STD),
    %% Generate signal
    ima = zeros(1,IMA_LEN);
    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+round((rand(1,1)-0.5)*IMA_LEN/2);
    ima(SIG_POS+[1:SIG_LEN]) = s; %% Insert the signal into the image array
    figure(i); subplot(3,1,1); plot(ima,'LineWidth',LN_WDTH);
    title('Image with signal and no noise.');
    
    
    %% Add noise to image
    ima = ima + NOISE_STD(i)*randn(1,IMA_LEN); %% This is white Gaussian noise
    figure(i); subplot(3,1,2); plot(ima,'LineWidth',LN_WDTH);
    title(['Image with white Gaussian noise of std dev \sigma=' num2str(NOISE_STD(i))]);
    
    %%
    %% Match-filter the image and calculate an appropriate threshold
    %%
    k = fliplr(s);           %% Set the kernel to be the matched filter
    for j=1:KSMOOTH, k = conv2(k,[1,2,1]/4); end;
    %% 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(i); 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');
end;