
%N=lunghezza del predittore
%alfa = eventuale ritardo
%b1 = i beta per la trasformazione dei tempi
%b2 = i polinomi per il bias space poli negativi
%b3 = i polinomi per il bias space poli positivi
%l = le varianze delle r.i.
%A = le relazioni inp/out definite da uscite e ingressi passati
%o2 = varianza innovazione
%Nb = numero di blocchi in cui partizionare i dati per scomporre Ey

function [C,Ey,Ey2]=Eybuild(N,alfa,b1,b2,b3,l,A,o2,BL,Nb)
warning off

m=size(A,3);

%%%%%%%%%%%%%%%%%%%%%%%%%lv%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if size(l,2)==1;
    lv=l*ones(1,m);
else
    lv=l;
end

%%%%%%%%%%%%%%%%%%%%%%%%%C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if size(b1,2)==1;
    beta=b1;
    n=N-alfa;
    t=abs(beta).^[1:n];%exp(-beta*[1:n])
    t=fliplr(t);
    for i=1:n;
        B(i,i:n)=[(t(i).^2)/2].*(t(i:n)-t(i)/3);
        B(i:n,i)=B(i,i:n)';
    end
    B=rot90(B);B=rot90(B);
    BB=(1e-6)*eye(alfa);
    for i=1:m;
        C(:,:,i)=[BB zeros(alfa,N-alfa);zeros(n,alfa) B];
    end
else
    b1v=b1;
    for i=1:m;
        beta=b1v(i);
        n=N-alfa;
        t=abs(beta).^[1:n];%exp(-beta*[1:n])
        t=fliplr(t);
        for i=1:n;
            B(i,i:n)=[(t(i).^2)/2].*(t(i:n)-t(i)/3);
            B(i:n,i)=B(i,i:n)';
        end
        B=rot90(B);B=rot90(B);
        BB=(1e-6)*eye(alfa);
        C(:,:,i)=[BB zeros(alfa,N-alfa);zeros(n,alfa) B];
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%B2 e B3 (bias spaces)%%%%%%%%%%%%%%%%
if ~isempty(b2) & size(b2,1)==1;%un solo bias space con poli negativi
    a2=dimpulse([1 zeros(1,max(size(b2)))],[1 b2],n);
    for i=1:m;
    B2(:,:,i)=toeplitz(a2,zeros(1,n));
    end
elseif ~isempty(b2) & size(b2,1)>1;
    for i=1:m
        a2=dimpulse([1 zeros(1,max(size(b2(i,:))))],[1 b2(i,:)],n);
        B2(:,:,i)=toeplitz(a2,zeros(1,n));
    end
end

if ~isempty(b3) & size(b3,1)==1;
    a3=dimpulse([1 zeros(1,max(size(b3)))],[1 b3],n);
    for i=1:m;
        B3(:,:,i)=toeplitz(a3,zeros(1,n));
    end
elseif ~isempty(b3) & size(b3,1)>1;
    for i=1:m
        a3=dimpulse([1 zeros(1,max(size(b3(i,:))))],[1 b3(i,:)],n);
        B3(:,:,i)=toeplitz(a3,zeros(1,n));
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%Ey%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if Nb==1;%una sola Ey da costruire
Ey=zeros(size(A,1),size(A,1));
for i=1:m;
    if ~isempty(b2) & ~isempty(b3);
        C(:,:,i)=lv(i)*B2(:,:,i)*B3(:,:,i)*C(:,:,i)*B3(:,:,i)'*B2(:,:,i)';
    elseif isempty(b2) & ~isempty(b3)
        C(:,:,i)=lv(i)*B3(:,:,i)*C(:,:,i)*B3(:,:,i)';
    elseif ~isempty(b2) & isempty(b3)
        C(:,:,i)=lv(i)*B2(:,:,i)*C(:,:,i)*B2(:,:,i)';
    elseif isempty(b2) & isempty(b3)
        C(:,:,i)=lv(i)*C(:,:,i);
    end
    Ey=Ey+A(:,:,i)*C(:,:,i)*A(:,:,i)';
end
Ey=Ey+o2*eye(size(A,1));
Ey2=[];
else%Ey a blocchi da restituire
    for kk=1:Nb;
        if kk<Nb;%primi blocchi
        I=[(1+(kk-1)*BL):(kk*BL)];
        Ey(:,:,kk)=zeros(size(I,2),size(I,2));
        for i=1:m;
            if ~isempty(b2) & ~isempty(b3)
                Ey(:,:,kk)=Ey(:,:,kk)+lv(i)*A(I,:,i)*B2(:,:,i)*B3(:,:,i)*C(:,:,i)*B3(:,:,i)'*B2(:,:,i)'*A(I,:,i)';
            elseif isempty(b2) & ~isempty(b3)
                Ey(:,:,kk)=Ey(:,:,kk)+lv(i)*A(I,:,i)*B3(:,:,i)*C(:,:,i)*B3(:,:,i)'*A(I,:,i)';
            elseif ~isempty(b2) & isempty(b3)
                Ey(:,:,kk)=Ey(:,:,kk)+lv(i)*A(I,:,i)*B2(:,:,i)*C(:,:,i)*B2(:,:,i)'*A(I,:,i)';
            elseif isempty(b2) & isempty(b3)
                Ey(:,:,kk)=Ey(:,:,kk)+lv(i)*A(I,:,i)*C(:,:,i)*A(I,:,i)';
            end
        end
        Ey(:,:,kk)=Ey(:,:,kk)+o2*eye(size(Ey(:,:,kk),1));
        else%ultimo blocco
            I=[(1+(kk-1)*BL):max(size(A,1))];
            Ey2=zeros(size(I,2),size(I,2));
            for i=1:m;
                if ~isempty(b2) & ~isempty(b3)
                    Ey2=Ey2+lv(i)*A(I,:,i)*B2(:,:,i)*B3(:,:,i)*C(:,:,i)*B3(:,:,i)'*B2(:,:,i)'*A(I,:,i)';
                elseif isempty(b2) & ~isempty(b3)
                    Ey2=Ey2+lv(i)*A(I,:,i)*B3(:,:,i)*C(:,:,i)*B3(:,:,i)'*A(I,:,i)';
                elseif ~isempty(b2) & isempty(b3)
                    Ey2=Ey2+lv(i)*A(I,:,i)*B2(:,:,i)*C(:,:,i)*B2(:,:,i)'*A(I,:,i)';
                elseif isempty(b2) & isempty(b3)
                    Ey2=Ey2+lv(i)*A(I,:,i)*C(:,:,i)*A(I,:,i)';
                end
            end
            Ey2=Ey2+o2*eye(size(Ey2,1));
        end
    end
end
    
