function [invLtpOutputBlock, invLtpOutputPast] = ...
    iltpfilt( invLtpInputBlock, subblockLength,ltpDelays, ltpGains, invLtpOutputPast )
% function [invLtpOutputBlock, invLtpOutputPast] = 
%     iltpfilt( invLtpInputBlock, subblockLength,ltpDelays, ltpGains, invLtpOutputPast )
%
% Open-Loop Inverse-LTP Filter
%
% invLtpInputBlock:   input vector
% subblockLength:     invLtpInputBlock is devided in subblocks of size subblockLength
% ltpDelays:          vector containing ltp delays (size: number of subblocks in invLtpInputBlock)
% ltpGains:           vector containing ltp gainss (size: number of subblocks in invLtpInputBlock)
% invLtpOutputPast:   inverse-ltp filter states before inverse-ltp filtering
% 
% invLtpOutputBlock:  inverse-ltp-filtered version of invLtpInputBlock
% invLtpOutputPast:   inverse-ltp filter states after inverse-ltp filtering
%  
% Author:  Markus Hauenstein
% Date:    03.01.2002 
% Contact: www.markus-hauenstein.de

blockLength = length(invLtpInputBlock);
maxDelay = length(invLtpOutputPast);
numberOfSubblocks = blockLength / subblockLength;

if numberOfSubblocks ~= floor(blockLength / subblockLength), 
  error( 'iltpfilt: blockLength and subblockLength do not fit' )
end

if numberOfSubblocks ~= length( ltpDelays )
  error( 'iltpfilt: number of delays and subblocks do not fit' )
end

if numberOfSubblocks ~= length( ltpGains )
  error( 'iltpfilt: number of gains and subblocks do not fit' )
end


invLtpOutputExt = [invLtpOutputPast; zeros(blockLength,1)];

startIndex = 1;
startIndexExt = maxDelay+1;
  
for j = 1:numberOfSubblocks,
  % (note that blockwise addition is not possible: delay might be < subblockLength )

  delayIndex = startIndexExt - ltpDelays(j);
  for k = 0:subblockLength-1,
    invLtpOutputExt(startIndexExt+k) = ... 
    invLtpInputBlock(startIndex+k) + ltpGains(j) * invLtpOutputExt(delayIndex+k); 
  end

  % set index for next subblock
  startIndex = startIndex + subblockLength;
  startIndexExt = startIndexExt + subblockLength;
end

invLtpOutputBlock = invLtpOutputExt( maxDelay+1 : maxDelay+blockLength );
invLtpOutputPast = invLtpOutputExt( blockLength+1 : maxDelay+blockLength );