function [cost_SID_DMST, key, p, shift, porm] = SID_DMST(ch,wl)

%Low-complexity filter design using shift-inclusive differential (SID) coefficients with DSMT (directed minimum spanining tree) solution
%
%Yongtao Wang, yongtaow@gmail.com; 
%
%Using SID, we have a greatly expanded design space, which is represented by a directed multiple graph
%The nodes (or vertices) of the graph is M=length(ch)+1. The M-th node is a virtual vertex.
%For all other nodes, each node represents a filter coefficient in ch, i.e. node j represents ch(j)
%In this graph, an edge from node M to another node, say j, is coefficient ch(j).
%There may be multiple edges directed from node i to j, (i,j not equal to M).
%Those edges represent SID coefficients ch(j)+(or -)2^L*ch(i). L is the amount of shift.
%Weight of an edge is defined as the number of adders required for computing that edge.
%A directed minimum spanning tree of the graph then yeilds a low-complexity implementation of the filter
%
%INPUTS:
%   ch: filter coefficients. For linear-phase filter, contains only the unique coefficients
%   wl: wordlength
%OUTPUTS:
%   cost_SID_DMST: number of adders required, i.e., total weight of the DMST
%   p: parent of each node from 1 to M. p(M)=0 since it is the root of DMST
%   key: weight of the edges in the DMST. key(i) is weight of the edge directed from p(i) to node i
%   shift: the amount of shift if a SID coefficient is used, see NOTE for more explanation
%   porm: "Plus OR Minus" in forming the SID coeff., see NOTE for more explanation
%NOTE: for j-th coefficient ch(j), if p(j)=M, then in implementing ch(j)*x (x is the inptu of the filter),
%   original value of ch(j) is used, i.e, no SID is used. In this case, shift(j)=0; porm(j)=0;
%   if p(j) is not M, then a SID coefficient ch(j)+porm(j)*2^shift(j)*ch(p(j)).
%   For example, let j=1, p(1)=2, shift(1)=-2, porm(1)=-1, then the SID coeff is ch(1)-2^(-2)*ch(2).
%   Hence, given p, shift and porm, we can easily determine the SID coeff used for the filter implementation
%KEY CAREABOUTS:
%1. edges directed from node i to j represent coefficients ch(j)+(or -)2^L*ch(i).Since the graph is directed, very 
%   important not to be confused over the direction of an edge and the coeff it represents
%2. elements of porm can be 1, -1 or 0. Then all the edges (i.e., coeff) directed from node p(j) to j
%   can be generally represented as ch(j)+porm(j)*2^shift(j)*ch(p(j)).
%3. key, p, shift and porm are also useful for debugging
%
%EXAMPLE:
% F = [0.3 .5]; A = [1 0]; wp=0.005; ws=0.005;
% DEV = [wp ws]; [n, fo, mo, w] = remezord(F,A,DEV);
% c = remez(n, fo, mo, w);
% wl = 16;
% cq=[];
% for i=1:length(c)
%     cq(i) = quantize(c(i), wl, 1);
% end
% ch = cq(1:ceil(length(c)/2));
%[cost_SID_DMST, key, p, shift, porm] = SID_DMST(ch,wl);

%construction of the directed multiple graph
M = length(ch)+1;
G = cell(M, M); GM={M,M};
shift=[]; porm=[];
%edges from cM to ci (i in [1, M-1]): original coefficients
for i=M
    for j=1:M-1
        edge = ch(j);  
        cbstr = CSD(edge, wl); %change CSD so that it outputs a vector consisting of 0, 1, and -1
        cost = sum(abs(cbstr));
        shift=0; %record the value of shift 
        porm=0; %plus or minus: 0 means that there is no differential coeff.
        if cost > 0
            cost = cost - 1;
        end
        G{j,i} = [edge; cost; shift; porm];
        ebpm = cbstr; %edge bit-pattern matrix
        GM{j,i} = ebpm;
    end
end
for i=1:M-1
    chb = CSD(ch(i), wl); ht = htdetect(chb);
    for j=1:M-1
        if i~=j
            edge = []; cost = []; shift=[]; porm=[];
            for L = 0 : ht(1)-1
                if abs(ch(j)-2^L*ch(i)) < 4/3
                    edge = [edge ch(j)-2^L*ch(i)];
                    porm=[porm -1];
                    shift = [shift L];
                end
                if abs(ch(j)+2^L*ch(i)) < 4/3
                    edge = [edge ch(j)+2^L*ch(i)];
                    porm=[porm 1];
                    shift = [shift L];
                end
            end
            for L = 1 : wl-ht(2)
                if abs(ch(j)-2^(-L)*ch(i)) < 4/3
                    edge = [edge ch(j)-2^(-L)*ch(i)];
                    porm=[porm -1];
                    shift = [shift -L];
                end
                if abs(ch(j)+2^(-L)*ch(i)) < 4/3
                    edge = [edge ch(j)+2^(-L)*ch(i)];
                    porm=[porm 1];
                    shift = [shift -L];
                end 
               
            end
            ebpm=[]; %edge bit pattern matrix
            for k = 1 : length(edge)
                cbstr = CSD(edge(k), wl);
                cost(k) = sum(abs(cbstr));
                ebpm(k,:) = cbstr;
            end
            GM{j,i} = ebpm;
            G{j,i} = [edge; cost;shift;porm]; 
        end
    end
end
%non-existent edges defined as edge=-100 and cost=1000
for i=1:M
    edge=[]; cost=[];
    edge = -100; cost = 1000; 
    shift=-100; porm=0;
    GM{i,i}=[];
    G{i,i} = [edge; cost;shift;porm];
end
for i=1:M-1
    edge=[]; cost=[];
    edge = -100; cost = 1000;
    shift=-100; porm=0;
    GM{M,i}=[];
    G{M,i} = [edge; cost;shift;porm];
end
%note: G_cost is a MxM matrix
%selected_edge_index: index of the selected edge 
G_cost=[]; selected_edge_index=[];
for i=1:M
    for j=1:M
        [Y, I] = min(G{i,j}(2,:));
        G_cost(i,j)=Y;
        selected_edge_index(i,j) = I(1);
    end
end
%DMST algorithm
[key, p] = DMST(G_cost);
%display('The cost with SID is:');
cost_SID_DMST = sum(key);

%Output shift and porm
shift=[]; porm=[];
for i=1:length(p)-1
    shift(i) = G{i,p(i)}(3,selected_edge_index(i,p(i)));
    porm(i) = G{i,p(i)}(4,selected_edge_index(i,p(i)));
end
