% SVD5c  Laguerre approximation
% For one interaction we will find the best fit.
% 

close all;
clear all;

%read data from database
[Resp,RPos0] = LeDadoIt(25);

DL = 47;
[a, b] = size(Resp');
%shapper for the signal that collects charge
Y=ones(a-4,1)-Resp(8,a+2-DL:a)';
Y=[zeros(4,1);Y];

%replace the contact that colected charge with the signal after going to the shapper
Resp(8,:)=Y';
Resp(:,1:DL)=Resp(:,a-DL+1:a);

%creates noise and add it to the input signal
noise = 0.00166 .* randn(size(Resp));
Respnoisy = Resp + noise;

%without noise
[Gains, yhat] = lagdecomp(Resp);
%with noise
[Gainsnoisy, yhatnoisy] = lagdecomp(Respnoisy);

% plots signal with and w/o noise
%for p=1:b
%    %without noise
%    figure, plot(Resp(p,:))';
%    hold;
%    plot(yhat(p,:),'r');
%    titletext = sprintf('signal %d of %d', p, b);
%    title(titletext);
%    %witnnoise
%    figure, plot(Respnoisy(p,:))';
%    hold;
%    plot(yhatnoisy(p,:),'g');
%    titletext = sprintf('signal with noise %d of %d', p, b);
%    title(titletext);
%end

%calculates the mse for both case on time
mse = diag((Resp-yhat)*(Resp-yhat)');
msewnoise = diag((Respnoisy-yhatnoisy)*(Respnoisy-yhatnoisy)');


for m=1:400
    %read data from database
    [Resp,RPos] = LeDadoIt(m);

    %shapper for the signal that collects charge
    Y=ones(a-4,1)-Resp(8,a+2-DL:a)';
    Y=[zeros(4,1);Y];

    %replace the contact that colected charge with the signal after going to the shapper
    Resp(8,:)=Y';
    Resp(:,1:DL)=Resp(:,a-DL+1:a);

    %without noise
    [Gains, yhat] = lagdecomp(Resp);
    
    mselaguerre(m) = sum(diag((Gains-Gainsnoisy)*(Gains-Gainsnoisy)'));
    
    msetime(m) = sum(diag((Resp-Respnoisy)*(Resp-Respnoisy)'));

end


p = [];
p2 = [];
for j=1:50
    %find minumum on both matrixes
    [x] = find(mselaguerre<=min(min(mselaguerre)))
    [x2] = find(msetime<=min(min(msetime)))
    
    mselaguerre(x) = 100;
    msetime(x2) = 100;
    
    %Do it for the segment that collected charge
    [Resp,RPos] = LeDadoIt(x);
    p = [p;RPos];
    d(j) = (RPos0-RPos)'*(RPos0-RPos);
    %Do it for all segments
    [Resp,RPos] = LeDadoIt(x2);
    p2 = [p2;RPos];
    d2(j) = (RPos0-RPos)'*(RPos0-RPos);
end
[p p2]




%sort mse for all elements
sortdmselaguerre = sort(mselaguerre);
figure, plot(sortdmselaguerre(1:5));
sortdmselaguerre(1:5)

%find minimum error and where it happend
noisemin = min(mselaguerre)
find(mselaguerre<=noisemin)
figure,plot(d(1:10))

%sort mse for all elements in time
sortdmsetime = sort(msetime);
figure, plot(sortdmsetime(1:5));
sortdmsetime(1:5)

figure, plot(10*sortdmselaguerre(1:12),'ko-');
hold;
plot(sortdmsetime(1:12),'kx-');
title('Distance of the points in order by mse from the original point calculated by time (black) and Laguerre methods (blue)')


%find minimum error and where it happend in time
noisemintime = min(msetime)
find(msetime <= noisemintime)
figure, plot(d(1:12),'ko-');
hold;
plot(d2(1:12),'kx-');
title('Distance of the points in order by mse from the original point calculated by time (black) and Laguerre methods (blue)')