function [S1,S2,Sig] = Projection2_(Sp_binary,STC_Matrix,Input_,p)
% Sp_Time are the spiking time represented by samples (not m-sec) 
L_v = size(STC_Matrix,1);
LL_V = 10*(L_v-1)+1; length_sta = LL_V;
Trial_num = size(Sp_binary,2);
LL = zeros(Trial_num,1);
for i =1:Trial_num
    spike_indx = find(Sp_binary(:,i)>0);
    maxindx = spike_indx(end);
    LL(i) = length(find(spike_indx>length_sta & spike_indx<maxindx-length_sta));
end
LL_ = [0;LL];
Sig = zeros(length_sta,sum(LL));
for i =1:Trial_num
    ind_=[];
    spike_indx = find(Sp_binary(:,i)>0);
    maxindx = spike_indx(end);
    ind_ = find(spike_indx>length_sta & spike_indx<maxindx-length_sta);
    Sig_ = zeros(length_sta,length(ind_));
    for j = 1:length(ind_)
        Sig_(:,j) = Input_(spike_indx(ind_(j))-length_sta+1:spike_indx(ind_(j)));
    end
    ind1 = sum(LL_(1:i)) + 1; ind2 = sum(LL_(1:i+1));
    Sig (:,ind1:ind2) = Sig_;
end
    
% [F_T, indx_T] = spike_times(Sp_binary,p);
% SP_Time = reshape(F_T,size(F_T,1)*size(F_T,2),1);
% Trial_num = size(Sp_binary,2);
%% Compute the eigenvectors and eigenvalues
% [u,s,v] = svd(STC_Matrix);
% %s_d = diag(sqrt(s)/sum(sqrt(diag(s))));%sqrt
% s_d = [1;1];%s/sum(diag(s));
% v1 = u(:,1);%s_d(1)*(u(:,1)-u(1,1));
% v2 = u(:,2);%s_d(2)*(u(:,2)-u(1,2));
% L_v = length(v1);
C = STC_Matrix;%cov(STC_Matrix');
[V,D] = eig(C);
d = diag(D)/sum(diag(D));
[q,p] = sort(d,'descend');
v1 = V(:,p(1));%sqrt(q(1))*
v2 = V(:,p(2));%sqrt(q(2))*
% v1 = sqrt(q(1))*(V(:,p(1))-V(1,p(1)));
% v2 = sqrt(q(2))*(V(:,p(2))-V(1,p(2)));
%% Resample to the originanl size
LL_V = 10*(L_v-1)+1;
L_F = LL_V-1;
uy1 = iddata([],v1,1);
ur1 = resample(uy1,LL_V,length(v1));
V1 = ur1.u;
uy2 = iddata([],v2,1);
ur2 = resample(uy2,LL_V,length(v2));
V2 = ur2.u;
%%
S1 = zeros(size(Sig,2),1);S2 = S1;
for j=1:length(S1)
S1(j) = Sig(:,j)'*V1/norm(Sig(:,j));
S2(j) = Sig(:,j)'*V2/norm(Sig(:,j));
end

% L_I = length(Input_);
% Inp_ = [zeros(L_F,1);Input_];%-mean(Input_)
% Sig = zeros(LL_V,length(find(SP_Time>0)));
% %length(SP_Time)
% j=0;
% for i = 1:round(length(SP_Time))
%     if SP_Time(i)~=0
%         j=j+1;
%     te = round(SP_Time(i)) + L_F;
%     Sig(:,j) = Inp_(te-L_F:te);
%     Sig(:,j) = Sig(:,j)/norm(Sig(:,j));
%     end
% end
% 
% s1 = Sig'*V1; %sqrt(s(1,1))*
% s2 = Sig'*V2; %sqrt(s(2,2))*
% s1 = Sig(1:L_V+1,:)'*V1(1:L_V+1);
% s2 = Sig(1:L_V+1,:)'*V2(1:L_V+1);