% pfftcnv2(x,b)
function y = pfftcnv2(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(1,Lx);
hp = zeros(P,C);
Hp = zeros(P,C);
xk = zeros(1,C);
Xk = zeros(1,C);
Xkps = zeros(S*P,C);
yk = zeros(1,C);
Y  = zeros(1,C);
xzp = zeros(1,Lx+C-L);
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,
  SLp = (p-1)*Np;
  hp(p,1:Np) = h(SLp+1:SLp+Np);
end;

% hp(k) => Hp(k)
for p=1:P,
  Hp(p,1:C) = fft(hp(p,1:C),C);
end;

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

tic;

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

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

end;
toc;

SUBPLOT(2,1,1),
  plot(1:N,h);
  n = [1:N];
  XT = [Np*n];
  set(gca,'XTick',XT);
  grid;
  PlotTitle = sprintf('Filterkoeffizienten h(n), N=%d, P=%d, S=%d, NP=%d',N,P,S,Np);
  title(PlotTitle);
  xlabel('n');
  ylabel('h');

SUBPLOT(2,1,2), 
  plot(1:K*L,real(y(1:K*L)));
  n = [1:N];
  XT = [L*n];
  set(gca,'XTick',XT);
  grid;
  PlotTitle = sprintf('Filterausgang y(n), K=%d, L=%d, C=%d',K,L,C);
  title(PlotTitle);
  xlabel('n');
  ylabel('y');