function y = ol_celp( x )
% function y = ol_celp( x )
%
% Simulation of CELP Codec with Open-Loop Pitch Predictor 
%
% x: codec input signal (sampled with 8 kHz, 16-bit quantized, i.e range -32768 ... 32767)
% y: codec output signal (sampled with 8 kHz, 16-bit quantized, i.e range -32768 ... 32767)
%   
% Author:  Markus Hauenstein
% Date:    16.01.2002 
% Contact: www.markus-hauenstein.de
%
% This program is copyrighted.

test = 0;

fprintf(1,'\nSimulation of CELP Codec with Open-Loop Pitch Predictor')
fprintf(1,'\nInitializing')


%
% init parameters identical for both coder and decoder which influence the bitrate
%

% ltp (number of samples)
numberOfSubblocks = 4;
subblockLength = 40;
minDelay = 20;
maxDelay = 147;
ltpGainsQuantBits = 6;

% lpc: bits for quantization of predictor coefficients
predQuantBits = 4*[1 1 1 1 1 1 1 1 1 1];

% size of stochastic gain codebook
stochGainCBSize = 64;

% size of stochastic Vec codebook
stochVecCBSize = 1024;


%
% init parameters identical for both coder and decoder which do NOT influence the bitrate
%

blockLength = numberOfSubblocks * subblockLength;

% ltp
ltpMaxGainQ = 1;

% lpc
lpcOrder = length(predQuantBits);
maxLarQ = 2.5;

% create stochastic gain codebook
minStochGain = -512;
maxStochGain = 512;
stochGainCBmu = 64;
stochGainCB = crmucb(minStochGain, maxStochGain, stochGainCBmu, stochGainCBSize);

% create normalized stochastic Vec codebook
stochVecCBAmp = 16;
randn('seed',0);
stochVecCB = randn( blockLength, stochVecCBSize );
for i=1:stochVecCBSize,
  oneDivRms  = 1/sqrt((stochVecCB(:,i)' * stochVecCB(:,i)) / blockLength);
  stochVecCB(:,i) = oneDivRms * stochVecCB(:,i);
end
stochVecCB = stochVecCBAmp * stochVecCB;

% calculate and display required bitrate
sampleFrequency = 8000;
lpcBitRate = sum( predQuantBits ) * sampleFrequency / blockLength;
ltpDelaysBitRate = ceil( log(maxDelay-minDelay)/log(2) ) * sampleFrequency / subblockLength;
ltpGainsBitRate = ltpGainsQuantBits * sampleFrequency / subblockLength;
stochVecBitRate = ceil( log(stochVecCBSize)/log(2) ) * sampleFrequency / blockLength;
stochGainBitRate = ceil( log(stochGainCBSize)/log(2) ) * sampleFrequency / blockLength;
bitRate = lpcBitRate + ltpDelaysBitRate + ltpGainsBitRate + stochVecBitRate + stochGainBitRate;
fprintf(1, '\nBitrate: %d bps', bitRate)


% format Input signal
x=x(:); % make column Vector
inputLength = length(x);
rest = inputLength - blockLength * floor(inputLength/blockLength);
if  rest == 0,
   numberOfBlocks = inputLength / blockLength;
else
   numberOfBlocks = ceil( inputLength / blockLength );
   x = [x; zeros( blockLength-rest, 1 )];
end


% initialize data storage for transmitted parameters 
aChannel = zeros(lpcOrder, numberOfBlocks);
ltpGainsChannel = zeros(numberOfSubblocks, numberOfBlocks);
ltpDelaysChannel = zeros(numberOfSubblocks, numberOfBlocks);
stochVecChannel = zeros(1, numberOfBlocks);
stochGainChannel = zeros(1, numberOfBlocks);


%
% THIS IS CODER 
%

fprintf(1,'\nCoding ')

% initialize coder variables

% ltp 
zLpc = zeros(lpcOrder,1);
ltpInputPast = zeros(maxDelay,1);
ltpGains = zeros(numberOfSubblocks,1);
ltpGainsQ = zeros(numberOfSubblocks,1);
ltpDelays = zeros(numberOfSubblocks,1);

% inverse ltp
invLtpOutputPast = zeros(maxDelay,1);

% inverse lpc 
zInvLpc = zeros(lpcOrder,1);

% Vec quantization of LPC/LTP residual - helper variables for Vec quantization
synFilterOutSignals = zeros( blockLength, stochVecCBSize );
zeroOutput = zeros(blockLength,1);
zeroInput = zeros(blockLength,1);

% coder main loop
blockIndex = 1 : blockLength;

for i=1:numberOfBlocks,
  fprintf(1,'.')

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % get current block
  xBlock = x( blockIndex );
  blockIndex = blockIndex + blockLength;


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % lpc (inverse vocal tract filter)
  lpcInputBlock = xBlock;

  % calculate linear prediction coefficients
  a = linpred( lpcInputBlock, lpcOrder );

  % simulate log-area-ratio quantization of LPC parameters
  aQ = larpredq( a, maxLarQ, predQuantBits );

  % lpc filter 
  [lpcOutputBlock, zLpc] = filter([1;-aQ], 1, lpcInputBlock, zLpc);


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % ltp (inverse pitch filter)
  ltpInputBlock = lpcOutputBlock;

  % determine LTP parameters
  [ltpDelays, ltpGains] = olltpdg( ltpInputBlock, subblockLength, minDelay, ltpInputPast );

  % simulate quantization of LTP parameters   
  ltpGainsQ = muquant( ltpGains, ltpMaxGainQ, ltpGainsQuantBits );
  
  % calculate ltp filter states for next looping
  [dummy, ltpInputPast] = ...
      ltpfilt( ltpInputBlock, subblockLength, ltpDelays, ltpGainsQ, ltpInputPast);


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Analysis-by-Synthesis: 
  % take LPC filter states from last looping into account and calculate a target vector
  zeroOutput = filter(1, [1;-aQ], zeroInput, zInvLpc);
  target = xBlock - zeroOutput; 


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Analysis-by-Synthesis: 
  % Synthesize all possible speech vectors by inverse ltp & lpc filtering 
  % of all noise vectors in stochastic codebook
  for l=1:stochVecCBSize,
    % Inverse LTP
    invLtpInputBlock = stochVecCB(:,l); 
    [invLtpOutputBlock, dummy] = ... 
        iltpfilt( invLtpInputBlock, subblockLength, ltpDelays, ltpGainsQ, invLtpOutputPast );
    % Inverse LPC
    invLpcInputBlock = invLtpOutputBlock;
    synFilterOutSignals(:,l) = filter( 1, [1;-aQ], invLpcInputBlock );
  end


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Analysis-by-Synthesis: 
  % Vector quantization of target vector   
  [bestStochGainIndex, bestStochVecIndex] = gsvecq( target, stochGainCB, synFilterOutSignals );


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Prepare next block by calculating filter states

  % calculate inverse ltp filter (pitch filter) states for next looping -
  % this corresponds to a decoder!

  invLtpInputBlock = stochVecCB(:,bestStochVecIndex); 
  [invLtpOutputBlock,invLtpOutputPast ] = ... 
      iltpfilt( invLtpInputBlock, subblockLength, ltpDelays, ltpGainsQ, invLtpOutputPast );

  % calculate inverse lpc filter (vocal tract filter) states
  invLpcInputBlock = stochGainCB(bestStochGainIndex) * invLtpOutputBlock;
  [dummy, zInvLpc] = filter(1, [1;-aQ], invLpcInputBlock, zInvLpc);

if test == 1,
  clg
  plot(xBlock,'g')
  hold
  plot(dummy,'r')
  pause
end


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % simulate transmission of parameters via channel to decoder
  aChannel(:,i) = aQ;
  ltpGainsChannel(:,i) = ltpGainsQ;
  ltpDelaysChannel(:,i) = ltpDelays;
  stochGainChannel(i) = bestStochGainIndex; 
  stochVecChannel(i) = bestStochVecIndex;

end  % coder main loop


%
% THIS IS DECODER 
%

fprintf(1,'\nDecoding')

% initialize decoder variables
y = zeros(numberOfBlocks*blockLength,1);

% inverse ltp
invLtpOutputPast = zeros(maxDelay, 1); 

% inverse lpc
zInvLpc = zeros(lpcOrder,1);

% decoder main loop
blockIndex = 1 : blockLength;

for i=1:numberOfBlocks,
  fprintf(1,'.')

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % inverse ltp (pitch reconstruction)

  %invLtpInputBlock = ltpGainBlockChannel(i) * stochVecCB(:,stochVecChannel(i)); 
  invLtpInputBlock = stochVecCB(:,stochVecChannel(i)); 

  [invLtpOutputBlock, invLtpOutputPast] = ... 
      iltpfilt( invLtpInputBlock, subblockLength, ltpDelaysChannel(:,i), ...
                ltpGainsChannel(:,i), invLtpOutputPast );


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % inverse lpc (vocal tract model)
  a = aChannel(:,i);
  invLpcInputBlock = invLtpOutputBlock;
  [invLpcOutputBlock, zInvLpc] = filter(1, [1;-a], invLpcInputBlock, zInvLpc);

  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % scaling
  yBlock = stochGainCB(stochGainChannel(i)) * invLpcOutputBlock; 


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % construct output signal
  y( blockIndex ) = yBlock;

  blockIndex = blockIndex + blockLength;
end  % decoder main loop


% format Output signal
y = y(1:inputLength);

if test == 1,
  clg
  plot(x,'g')
  hold
  plot(y,'r')
end

fprintf(1,'\nReady\n')