% pflmss(mu,N,C,P,S,L,x,d,w_start)
function [ee,w,dd,calc,PX]=pflmss(mu,N,C,P,S,L,x,d,w_start)

Lx = length(x);
Np = S*L;
K = Lx/L;

% Initialisierungen

%% FLMS Parameter
gamma = 0.6;	             % Vergessensfaktor

ee=zeros(Lx,1);     % Fehlervektor
dd=zeros(Lx,1);     % Filterausgang, Schaetzung von d
wp = zeros(C,P);	% Gewichtsvektor im Frequenzbereich
w  = zeros(P*S*L,1);        % Gewichtsvektor im Zeitbereich
PX = ones(C,1);	% Schaetzung des  Leistungsdichtespektrums

WSp = zeros(C,P);
X = zeros(C,S*P);
U = zeros(C,S*P);
ys = zeros(C,1);
YS  = zeros(C,1);
xzp = zeros(Lx+S*L,1);
pidTbl = [2:P*S];
pidTbl(P*S) = 1;
pStbl = [0:P*S-1];
pStbl(1) = P*S;
pjiTbl = [2:P];
pjiTbl(P) = 1;

xzp(C-L+1:Lx+C-L) = x;

tic;     % Stopuhr laeuft

pid = 1;
pji = 1;

flops(0);

for k=1:K,
  
  kL = (k-1)*L;
  
  X(1:C,pid) = fft(xzp(kL+1:kL+C),C);
    
  pS = pid;
  YS = WSp(1:C,1).*X(1:C,pid);
  for p=2:P,
    for i = 1:S
      pS = pStbl(pS);
    end;
    YS   = YS + WSp(1:C,p).*X(1:C,pS);
  end;

  ys = real(ifft(YS,C));
  dd(kL+1:kL+L) = ys(C-L+1:C);
  
  ee(kL+1:kL+L) = d(kL+1:kL+L) - ys(C-L+1:C);
  E  = fft([zeros(C-L,1); ee(kL+1:kL+L)]);
  PX = abs((1-gamma)*conj(X(1:C,pid)).*X(1:C,pid) + gamma*PX);
  
  U(1:C,pid) = mu ./ (PX+0.001);
  
  % Update
  pS = pid;
  for p=1:P,
    WSp(1:C,p) = WSp(1:C,p) + U(1:C,pS) .* conj(X(1:C,pS)) .* E;
    for i = 1:S
      pS = pStbl(pS);
    end;
  end;

% Teuer: Projektion jeden p-ten Teilfilters wp pro k-ter Iteration
%  for p=1:P,
%    wp(1:C,p)  = real(ifft(WSp(1:C,p)));
%    WSp(1:C,p) = fft(wp(1:Np,p),C);
%  end;

% Billig: Projektion eines Teilfilters alternierend pro k-ter Iteration
  wp(1:C,pji)  = real(ifft(WSp(1:C,pji)));
  WSp(1:C,pji) = fft(wp(1:Np,pji),C);
  
  pji = pjiTbl(pji);
  pid = pidTbl(pid);

end;

flops

% Return estimated filter weights
for p=1:P,
  w((p-1)*Np+1:p*Np) = wp(1:Np,p);
end;

calc = toc;