clear
close all
N=4000; %data size
%v=[0:N-1];
fs=30;
v=0:1/fs:(4000-1)/fs;
t=v;
%s_i=100.*((rem(v,50)-13)/9)';
xt=[0:9 10:-1:1];
%s_i =100*repmat(xt,1,200);
s_i=100.*sin(2.*3*pi.*v); %+80*cos(2.*12*pi.*v)+120*cos(2.*20*pi.*v)+100.*((rem(v,5)-13)/9);
s_i=s_i';
s_in=s_i;%+50*randn(size(s_i));
figure
plot(v,s_in)
sr1=s_in;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Analisi Dinamica %%%%%%%%%%%%%%%%%%%
s1=signal(sr1);
figure
subplot(2,2,1), view(s1)
  grid
a1 = amutual(s1,128);
%laga1=firstmin(a1);
subplot(2,2,2), view(a1)
  grid
 laga1=input('Choose the time lag: first minimum of amutual corr.  ');
c1 = cao(s1,30,laga1,10,500);
subplot(2,2,3), view(c1)
  grid
%cmax1=firstmax(c1);
cmax1=input('Choose the embedding dimension  ');
e1r = embed(s1, 3, laga1);
e1r2=embed(s1, 2, laga1);
subplot(2,2,4),view(e1r)
figure
view(e1r2)
e1 = embed(s1, cmax1, laga1); %obviously we see nothing
st1=data(e1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Recurrence Plot 
srP=sr1(1:10:end);m = cmax1; tau = laga1;
trP=t(1:10:end)';
N = length(srP);
N2 = N - tau * (m - 1);
%The original data series has a length of 600, but the resulting phase space trajectory has the length 
%596. Now we create the phase space trajectory with 
%clear xe
for mi = 1:m;
    xe(:, mi) = srP([1:N2] + tau * (mi-1));
end
%We can accelerate the pair-wise test between each points on the trajectory with a fully vectorized 
%algorithm. For that we need to transfer the trajectory vector into two test vectors, whose 
%component-wise test will provide the pair-wise test of the trajectory vector: 
x1 = repmat(xe, N2, 1);
x2 = reshape(repmat(xe(:), 1, N2)', N2 * N2, m);
%Using these vectors we calculate the recurrence plot using the Euclidean norm without any loop
S = sqrt(sum( (x1 - x2) .^ 2, 2 ));
S = reshape(S, N2, N2);
figure
imagesc(trP(1:N2), trP(1:N2), S < 25)
axis square
colormap([1 1 1;1 0 0])
xlabel('Time (sec)'), ylabel('Time (sec)')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%L=length(s_in)/2;
%fs=32;
 pause
 L=cmax1;
I=[1:cmax L];
x1=s_in;
% Step1 : Build trajectory matrix
   N=length(x1); 
   if L>N/2;L=N-L;end
	K=N-L+1; 
   X=zeros(L,K);  
	for i=1:K
	  X(1:L,i)=x1(i:L+i-1); 
    end
    figure
   plot(X(1,:),X(2:cmax1,:),'-o');
% Step 2: SVD

   S=X*X'; 
	[U,autoval]=eig(S);
	[d,i]=sort(-diag(autoval));  
   d=-d;
   U=U(:,i);sev=sum(d); 
   figure
	loglog((d./sev)*100),hold on,loglog((d./sev)*100,'ro');
%    xlim([0.5 L+50]); ylim([0.001 200]);
	title('Singular Spectrum');xlabel('Eigenvalue Number');ylabel('Eigenvalue (% Norm of trajectory matrix retained)')
   grid
    V=(X')*U; 
   rc=U*V';

% Step 3: Grouping

  % I=input('Choose the agrupation of components to reconstruct the series in the form I=[i1,i2:ik,...,iL];  ');
   Vt=V';
   rca=U(:,I)*Vt(I,:);

% Step 4: Reconstruction

   y=zeros(N,1);  
   Lp=min(L,K);
   Kp=max(L,K);

   for k=0:Lp-2
     for m=1:k+1;
      y(k+1)=y(k+1)+(1/(k+1))*rca(m,k-m+2);
     end
   end

   for k=Lp-1:Kp-1
     for m=1:Lp;
      y(k+1)=y(k+1)+(1/(Lp))*rca(m,k-m+2);
     end
   end

   for k=Kp:N
      for m=k-Kp+2:N-Kp+1;
       y(k+1)=y(k+1)+(1/(N-k))*rca(m,k-m+2);
     end
   end

   figure;subplot(2,1,1);hold on;xlabel('Data poit');ylabel('Original and reconstructed series')
   plot(x1);grid on;plot(y,'r')
% size(x1)
% size(y)
   r=x1-y;
   subplot(2,1,2);plot(r,'g');xlabel('Data poit');ylabel('Residual series');grid on
   vr=(sum(d(I))/sev)*100;
   [Wefica1, ISRef1, Wsymm1,ISRsymm1, status1, icasig_E1]=efica_t(st1');
%[Wefica, ISRef, Wsymm,Werr, ISRsymm, status, icasig_E]=efica_t1(st');
assignin('base', 'myplot_signal', icasig_E1);
myplot
 fre1 = figure('Visible', 'on', 'Color', [0.753 0.753 0.753], 'Name', 'ISR matrix of EFICA', 'Resize', 'off', 'MenuBar', 'none', 'NumberTitle', 'off');
          showmatrixISR(ISRsymm1);