function h3d = diff3(h3,N,a1)
%DIFF3 Difference of triangular 3rd order kernel in vector form
%
% H3D = DIFF3(H3,N,A) returns the triangular kernel associated to
% the generalized transfer function 
% H3D(z1,z2,z3)=(1-A/z1)(1-A/z2)(1-A/z3)H3(z1,z2,z3)
%
% N is the memory length (in samples) of H3 and H3D, which are vectorized

%   CBVK Toolbox
%   JHMG, PMSB 2016


h3d = zeros(size(h3));

a2 = a1*a1;
a3 = a1*a2;

ind = 0;

ind = ind+1;
% (1,1,1)
h3d(ind) = h3(ind);

ind = ind+1;
% (1,1,2)
h3d(ind) =  - a1*h3( ind + -1 ) + (1/3)*h3(ind);
h3d(ind) = 3*h3d(ind);

for k3=1+2:N
    ind = ind+1;
    % (1,1,k3)
    h3d(ind) =  - a1*(1/3)*h3( ind + -1 ) + (1/3)*h3(ind);
    h3d(ind) = 3*h3d(ind);
    
end

ind = ind+1;
% (1,2,2)
h3d(ind) = a2*h3( ind + -N ) - a1*2*(1/3)*h3( ind + 1 - N ) + (1/3)*h3(ind);
h3d(ind) = 3*h3d(ind);

ind = ind+1;
% (1,2,3)
h3d(ind) = a2*(1/3)*h3( ind + -N ) - a1*(1/3)*h3( ind + -1 ) - a1*(1/3)*h3( ind + 1 - N ) + (1/6)*h3(ind);
h3d(ind) = 6*h3d(ind);

for k3=2+2:N
    ind = ind+1;
    % (1,2,k3)
    h3d(ind) = a2*(1/3)*h3( ind + -N ) - a1*(1/6)*h3( ind + -1 ) - a1*(1/3)*h3( ind + 1 - N ) + (1/6)*h3(ind);
    h3d(ind) = 6*h3d(ind);
    
end

for k2=1+2:N
    ind = ind+1;
    % (1,k2,k2)
    h3d(ind) = a2*(1/3)*h3( ind + k2 - N - 2 ) - a1*2*(1/6)*h3( ind + k2 - N - 1 ) + (1/3)*h3(ind);
    h3d(ind) = 3*h3d(ind);
    
    
    if (k2 <= N-1)
        % (1,k2,k2+1)
        ind = ind+1;
        h3d(ind) = a2*(1/6)*h3( ind + k2 - N - 2 ) - a1*(1/3)*h3( ind + -1 ) - a1*(1/6)*h3( ind + k2 - N - 1 ) + (1/6)*h3(ind);
        h3d(ind) = 6*h3d(ind);
        
    end
    for k3=k2+0+2:N
        ind = ind+1;
        % (1,k2,k3)
        h3d(ind) = a2*(1/6)*h3( ind + k2 - N - 2 ) - a1*(1/6)*h3( ind + -1 ) - a1*(1/6)*h3( ind + k2 - N - 1 ) + (1/6)*h3(ind);
        h3d(ind) = 6*h3d(ind);
        
    end
end

for k1=2:N
    ind = ind+1;
    % (k1,k1,k1)
    k2 = k1;
    h3d(ind) =  - a3*h3( ind + (3*k1)/2 - (5*N)/2 + k2 + N*k1 - N^2/2 - k1^2/2 - 3 ) + a2*3*(1/3)*h3( ind + (3*k1)/2 - (5*N)/2 + k2 + N*k1 - N^2/2 - k1^2/2 - 2 ) - a1*3*(1/3)*h3( ind + -((N - k1 + 1)*(N - k1 + 2))/2 ) + h3(ind);
    
    
    if (k1 <= N-1)
        % (k1,k1,k1+1)
        ind = ind+1;
        h3d(ind) =  - a3*(1/3)*h3( ind + (3*k1)/2 - (5*N)/2 + k2 + N*k1 - N^2/2 - k1^2/2 - 3 ) + a2*2*(1/3)*h3( ind + (3*k1)/2 - (3*N)/2 + N*k1 - N^2/2 - k1^2/2 - 2 ) - a1*h3( ind + -1 ) + a2*(1/3)*h3( ind + (3*k1)/2 - (5*N)/2 + k2 + N*k1 - N^2/2 - k1^2/2 - 2 ) - a1*2*(1/6)*h3( ind + -((N - k1 + 1)*(N - k1 + 2))/2 ) + (1/3)*h3(ind);
        h3d(ind) = 3*h3d(ind);
        
    end
    for k3=k1+0+2:N
        ind = ind+1;
        % (k1,k1,k3)
        h3d(ind) =  - a3*(1/3)*h3( ind + (3*k1)/2 - (5*N)/2 + k2 + N*k1 - N^2/2 - k1^2/2 - 3 ) + a2*2*(1/6)*h3( ind + (3*k1)/2 - (3*N)/2 + N*k1 - N^2/2 - k1^2/2 - 2 ) - a1*(1/3)*h3( ind + -1 ) + a2*(1/3)*h3( ind + (3*k1)/2 - (5*N)/2 + k2 + N*k1 - N^2/2 - k1^2/2 - 2 ) - a1*2*(1/6)*h3( ind + -((N - k1 + 1)*(N - k1 + 2))/2 ) + (1/3)*h3(ind);
        h3d(ind) = 3*h3d(ind);
        
    end
    
    if (k1 <= N-1 && k1 <= N-1)
        % (k1,k1+1,k1+1)
        k2 = k1+1;
        ind = ind+1;
        h3d(ind) =  - a3*(1/3)*h3( ind + (3*k1)/2 - (5*N)/2 + k2 + N*k1 - N^2/2 - k1^2/2 - 3 ) + a2*h3( ind + k2 - N - 2 ) + a2*2*(1/6)*h3( ind + (3*k1)/2 - (5*N)/2 + k2 + N*k1 - N^2/2 - k1^2/2 - 2 ) - a1*2*(1/3)*h3( ind + k2 - N - 1 ) - a1*(1/3)*h3( ind + -((N - k1 + 1)*(N - k1 + 2))/2 ) + (1/3)*h3(ind);
        h3d(ind) = 3*h3d(ind);
        
    end
    
    if (k1 <= N-1 && k1 <= N-2)
        % (k1,k1+1,k1+2)
        ind = ind+1;
        h3d(ind) =  - a3*(1/6)*h3( ind + (3*k1)/2 - (5*N)/2 + k2 + N*k1 - N^2/2 - k1^2/2 - 3 ) + a2*(1/3)*h3( ind + k2 - N - 2 ) + a2*(1/3)*h3( ind + (3*k1)/2 - (3*N)/2 + N*k1 - N^2/2 - k1^2/2 - 2 ) - a1*(1/3)*h3( ind + -1 ) + a2*(1/6)*h3( ind + (3*k1)/2 - (5*N)/2 + k2 + N*k1 - N^2/2 - k1^2/2 - 2 ) - a1*(1/3)*h3( ind + k2 - N - 1 ) - a1*(1/6)*h3( ind + -((N - k1 + 1)*(N - k1 + 2))/2 ) + (1/6)*h3(ind);
        h3d(ind) = 6*h3d(ind);
        
    end
    for k3=k1+1+2:N
        
        if (k1 <= N-1)
            % (k1,k1+1,k3)
            ind = ind+1;
            h3d(ind) =  - a3*(1/6)*h3( ind + (3*k1)/2 - (5*N)/2 + k2 + N*k1 - N^2/2 - k1^2/2 - 3 ) + a2*(1/3)*h3( ind + k2 - N - 2 ) + a2*(1/6)*h3( ind + (3*k1)/2 - (3*N)/2 + N*k1 - N^2/2 - k1^2/2 - 2 ) - a1*(1/6)*h3( ind + -1 ) + a2*(1/6)*h3( ind + (3*k1)/2 - (5*N)/2 + k2 + N*k1 - N^2/2 - k1^2/2 - 2 ) - a1*(1/3)*h3( ind + k2 - N - 1 ) - a1*(1/6)*h3( ind + -((N - k1 + 1)*(N - k1 + 2))/2 ) + (1/6)*h3(ind);
            h3d(ind) = 6*h3d(ind);
            
        end
    end
    for k2=k1+0+2:N
        ind = ind+1;
        % (k1,k2,k2)
        h3d(ind) =  - a3*(1/3)*h3( ind + (3*k1)/2 - (5*N)/2 + k2 + N*k1 - N^2/2 - k1^2/2 - 3 ) + a2*(1/3)*h3( ind + k2 - N - 2 ) + a2*2*(1/6)*h3( ind + (3*k1)/2 - (5*N)/2 + k2 + N*k1 - N^2/2 - k1^2/2 - 2 ) - a1*2*(1/6)*h3( ind + k2 - N - 1 ) - a1*(1/3)*h3( ind + -((N - k1 + 1)*(N - k1 + 2))/2 ) + (1/3)*h3(ind);
        h3d(ind) = 3*h3d(ind);
        
        
        if (k2 <= N-1)
            % (k1,k2,k2+1)
            ind = ind+1;
            h3d(ind) =  - a3*(1/6)*h3( ind + (3*k1)/2 - (5*N)/2 + k2 + N*k1 - N^2/2 - k1^2/2 - 3 ) + a2*(1/6)*h3( ind + k2 - N - 2 ) + a2*(1/3)*h3( ind + (3*k1)/2 - (3*N)/2 + N*k1 - N^2/2 - k1^2/2 - 2 ) - a1*(1/3)*h3( ind + -1 ) + a2*(1/6)*h3( ind + (3*k1)/2 - (5*N)/2 + k2 + N*k1 - N^2/2 - k1^2/2 - 2 ) - a1*(1/6)*h3( ind + k2 - N - 1 ) - a1*(1/6)*h3( ind + -((N - k1 + 1)*(N - k1 + 2))/2 ) + (1/6)*h3(ind);
            h3d(ind) = 6*h3d(ind);
            
        end
        for k3=k2+0+2:N
            ind = ind+1;
            % (k1,k2,k3)
            h3d(ind) =  - a3*(1/6)*h3( ind + (3*k1)/2 - (5*N)/2 + k2 + N*k1 - N^2/2 - k1^2/2 - 3 ) + a2*(1/6)*h3( ind + k2 - N - 2 ) + a2*(1/6)*h3( ind + (3*k1)/2 - (3*N)/2 + N*k1 - N^2/2 - k1^2/2 - 2 ) - a1*(1/6)*h3( ind + -1 ) + a2*(1/6)*h3( ind + (3*k1)/2 - (5*N)/2 + k2 + N*k1 - N^2/2 - k1^2/2 - 2 ) - a1*(1/6)*h3( ind + k2 - N - 1 ) - a1*(1/6)*h3( ind + -((N - k1 + 1)*(N - k1 + 2))/2 ) + (1/6)*h3(ind);
            h3d(ind) = 6*h3d(ind);
            
        end
    end
end

