% pfftfilt(x,b)
function y = pfftfilt(x,h,P,S,L)

N = length(h);
Lx = length(x);
C = 2^ceil(log(L+S*L-1)/log(2));
Np = S*L;
K = Lx/L;

y = zeros(Lx,1);
Hp = zeros(C,P);
Xkps = zeros(C,S*P);
yk = zeros(C,1);
Y  = zeros(C,1);
xzp = zeros(Lx+C-L,1);
pidTbl = [2:P*S];
pidTbl(P*S) = 1;
pStbl = [0:P*S-1];
pStbl(1) = 10;

pid = 1;

% h(k) => Hp(k)
for p=1:P,
  Hp(1:C,p) = fft([h((p-1)*Np+1:p*Np);zeros(C-Np,1)],C);
end;

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

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

  yk = ifft(Y,C);
 
  y(kL+1:kL+L) = yk(C-L+1:C);
 
  pid = pidTbl(pid);

end;