% SVD5  Laguerre approximation
% Question: finding the MSE in time is the same as finding the
% MSE for the adjustable parameters????
% uses the samples of the segment that collects the charge,
% 

clear
hold off
format compact
T=.1; p = .5
K = sqrt((1-p^2)*T);
DL=47;
n = [0 K];               %find the outputs of the laguerre polynomial
d = [1 -1*p];
u = dimpulse(n,d,DL); 
nt = [-p 1];
dt = [1 -p];

%Creates the Laguere base with par waveforms and the size set by DL
par=10
for i=1:par-2
 u = [u dlsim(nt,dt,u(:,i))]; 
end
u = [ ones(47,1)./4.7 u];
u(1,1) = 0;

%u=[[zeros(3,16)];u];
%u=[[.3;zeros(46,1)],u];
%u=[[0;.3;zeros(45,1)],u];
%u=[[0;0;.3;zeros(44,1)],u];

%loads all shapes that were previously saved on a workspace
load shapes2
%shapper for the signal that collects charge
%Y=ones(46,19)-shapes(:,50+2-DL:50)';
%Y=[zeros(1,19);Y];
Y = shapes(:,50+2-DL:50)';
Y=[zeros(1,19);Y];
I=size(Y);

%creates Y_hat
y_hat=zeros(I);

%calculates  the inver os this matrix
invutuut=inv(u'*u)*u';

for i=1:I(1,2)
  %t are the Laguere vector of weights
  t(:,i)=invutuut*Y(:,i);
  %temporal rep of the function when going back from Laguere decomposition
  y_hat(:,i)=u*t(:,i);
  %error in time
  mse(i)=[Y(:,i)-y_hat(:,i)]'*[Y(:,i)-y_hat(:,i)];
%figure(4),plot(1:DL,Y(:,i),'r',1:DL,y_hat(:,i),'g')
%pause(1)
end

%sum(mse)

%measured
%creates a waveform the is the linear sum of two wave forms from the base plus noise
randn('seed',0);
M=2*Y(:,1)+3*Y(:,11)+.05*randn(47,1)-.025;

%   find mse for time sequence
%mse=.1*ones(I(1,2),I(1,2));
%for i=1:I(1,2)
% for j=1:I(1,2)
%  if i~=j
%   g=inv([Y(:,i)';Y(:,j)']*[Y(:,i),Y(:,j)])*[Y(:,i)';Y(:,j)']*M;  
%   mse(i,j)=[M-g(1)*Y(:,i)-g(2)*Y(:,j)]'*[M-g(1)*Y(:,i)-g(2)*Y(:,j)];
%figure(1);plot(1:47,g(1)*Y(:,i)+g(2)*Y(:,j),'g',1:47,M,'r')  %plot on time
%[i,j,mse(i,j)]
%pause(1)
%  end
% end
%end

%[r1,c1]=find(mse<=(min(min(mse))))

%find the Laguerre parameters for the measured signal.
Mt=invutuut*M;
M_hat=u*Mt;
%mse=[M-M_hat]'*[M-M_hat]
figure(5);plot(1:DL,M,'r',1:DL,M_hat,'g')


%   find mse for Laguerre parameter sequence
mset=1*ones(I(1,2),I(1,2));
mse=.1*ones(I(1,2),I(1,2));
for i=1:I(1,2)
 for j=1:I(1,2)
  if i~=j
%tempo
   g=inv([Y(:,i)';Y(:,j)']*[Y(:,i),Y(:,j)])*[Y(:,i)';Y(:,j)']*M;  
   mse(i,j)=[M-g(1)*Y(:,i)-g(2)*Y(:,j)]'*[M-g(1)*Y(:,i)-g(2)*Y(:,j)];
figure(1);plot(1:47,g(1)*Y(:,i)+g(2)*Y(:,j),'g',1:47,M,'r')  %plot on time

%Laguere
   gt=inv([t(:,i)';t(:,j)']*[t(:,i),t(:,j)])*[t(:,i)';t(:,j)']*Mt;  
   mset(i,j)=[Mt-gt(1)*t(:,i)-gt(2)*t(:,j)]'*[Mt-gt(1)*t(:,i)-gt(2)*t(:,j)];
figure(2);plot(1:47,M,'w',1:47,g(1)*Y(:,i)+g(2)*Y(:,j),'y',1:47,u*(gt(1)*t(:,i)+gt(2)*t(:,j)),'g',1:47,M_hat,'r')  %plot on time
figure(3),plot(1:par,gt(1)*t(:,i)+gt(2)*t(:,j),'g',1:par,Mt,'r')  %plot on Laguere
[i,j,mse(i,j),.1*mset(i,j),g',gt']

%pause(1)
  end 
 end
end

[r,c]=find(mse<=(min(min(mse))))
[rt,ct]=find(mset<=(min(min(mset))))

sum(sum(mse))
sum(sum(.1*mset))

%[z,w]=min(mse);
%[zt,wt]=min(mset);

