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