function [M,DOF] = make_delta_blocking(nu,N,dibtype)
%% Function to obtain Delta-input-blocking matrix
% nu - number of inputs
% N - prediction horizon
% dibtype - for how many inputs is difference constant for exampe dibtype =
% [2 3] => [(u1,u2) (u2,u3,u4)] => u2 - u1 = C1, u3 - u2 = C2, u4 - u3 = C2
ndibt = size(dibtype,2); % number of delta input blocking type
neqM = 0;  % number of equations matrix M (just for now 0)

for s = 1 : ndibt

NN = dibtype(s);    % number of prediction horizon for small matrix MM 
neqAA = NN-1;  % number of equations matrix AA
neqMM = NN-2;  % number of equations matrix MM
MM = cell(neqMM,NN);    % size MM
AA = cell(neqAA,NN);    % size AA


%% matrix AA 
% matrix where we define eg. u2 - u1 = C, u3 - u2 = C.....
% zeros
for i = 1:neqAA
    for j = 1:NN
        AA{i,j} = zeros(nu);
    end
end

% - ones
for i = 1:neqAA
    for j = i
        AA{i,j} = -eye(nu);
    end
end

% ones
for i = 1:neqAA
    for j = i+1
        AA{i,j} = eye(nu);
    end
end

%% matrix MM
% make matrix MM where we from AA remove constant C

for i = 1:neqMM
    for j = 1:NN
        MM{i,j} = AA{i,j}-AA{neqAA,j};  % remove constant C
    end
end

bigM{s}=MM; % big M consist of small MM matrix using NN as prediction horizon
neqM =  neqM + size(bigM{s},1);  % number of rows in M

end

%% delta input blocking matrix M

M = cell(neqM,N);   % size M

% zeros in M
for i = 1:neqM
    for j = 1:N
        M{i,j} = zeros(nu);
    end
end

% put MMs into M
sizeBMprew = [0 0]; % "starting point"
for s = 1 : ndibt
   
    sizeBM = size(bigM{s});
    for i = sizeBMprew(1)+1 : sizeBMprew(1)+sizeBM(1)
        for j = sizeBMprew(2)+1:sizeBMprew(2)+(sizeBM(2))
            M{i,j} = bigM{s}{i-sizeBMprew(1),j-(sizeBMprew(2))};
        end
    end
    sizeBMprew = sizeBM+(sizeBMprew-[0 1]);
end

M = cell2mat(M);

% degrees of freedom
DOF = N - neqM;


