function [Qf,Rf,Qb,Rb,V,S,L]=dfe(B,A,M,N,rho,m);
%function [Qf,Rf,Qb,Rb,V]=dfe(B,A,M,N,rho,m);
%
% optimal dfe with mse performance given correct past decision
% sb950315

if nargin<6,m=0;end

% matrix P
nb=length(B)-1;
b1=[B zeros(1,m+1-nb)];b1=b1(1:m+1);
na=length(A)-1;
a1=[A zeros(1,m+1-na)];a1=a1(1:m+1);
aa=toeplitz(a1,[1 zeros(1,m)]);
p=inv(aa)*b1';p=p';
P=toeplitz(p,[p(1) zeros(1,m)]);

% matrix T
nm=length(M)-1;
m1=[M zeros(1,m+1-nm)];m1=m1(1:m+1);
nn=length(N)-1;
n1=[N zeros(1,m+1-nn)];n1=n1(1:m+1);
aa=toeplitz(n1,[1 zeros(1,m)]);
t=inv(aa)*m1';t=t';
T=toeplitz(t,[t(1) zeros(1,m)]);

% gamma-values and polynomials S,L,L1
r=[zeros(m,1);1];
gamma=inv(P'*P+rho*T'*T)*P'*r;
S=T*gamma;S=S';
L1=-P*gamma;L1=L1';
L=L1+r';

%polynomial Q
aus1=conv(conv(B,N),S);
aus2=conv(conv(A,M),L1);
laus1=length(aus1);
laus2=length(aus2);
aus1=[aus1 zeros(1,laus2-laus1)];
aus2=[aus2 zeros(1,laus1-laus2)];
aus=aus1+aus2;
laus=max(laus1,laus2);
aus=[aus zeros(1,m+2-laus)];
laus=length(aus);
Q=aus(m+2:laus);

% feedback and feedforward filters
Qb=Q;
Rb=conv(A,M);
Qf=conv(S,N);
Rf=M;

% mse
V=sum(L.^2+rho*S.^2);