function y = cl_celp( x )
% function y = cl_celp( x )
%
% Simulation of CELP Codec with Closed-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 Closed-Loop Pitch Predictor')
fprintf(1,'\nInitializing')

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

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

% size of adaptive gain codebook
adapGainCBSize = 64;

% 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 Vector codebook
stochVecCBSize = 1024;


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

blockLength = numberOfSubblocks * subblockLength;

% adaptive Vector codebook
adapVecCBSize = maxDelay-minDelay+1;

% create adaptive gain codebook
minAdapGain = -1024;
maxAdapGain = 1024;
adapGainCBmu = 64;
adapGainCB = crmucb(minAdapGain, maxAdapGain, adapGainCBmu, adapGainCBSize);

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

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

% create normalized stochastic vector codebook
stochVecCBSize = 1024;
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;
adapVecBitRate = ceil( log(adapVecCBSize)/log(2) ) * sampleFrequency / subblockLength;
adapGainBitRate = ceil( log(adapGainCBSize)/log(2) ) * sampleFrequency / subblockLength;
stochVecBitRate = ceil( log(stochVecCBSize)/log(2) ) * sampleFrequency / blockLength;
stochGainBitRate = ceil( log(stochGainCBSize)/log(2) ) * sampleFrequency / blockLength;
bitRate = lpcBitRate + adapVecBitRate + adapGainBitRate + 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);
adapGainChannel = zeros(numberOfSubblocks, numberOfBlocks);
adapVecChannel = zeros(numberOfSubblocks, numberOfBlocks);
stochGainChannel = zeros(1, numberOfBlocks);
stochVecChannel = zeros(1, numberOfBlocks);


%
% THIS IS CODER 
%

z = zeros(numberOfBlocks*blockLength,1);


fprintf(1,'\nCoding ')

% initialize coder variables

% ltp 
bestAdapGainIndices = zeros(numberOfSubblocks,1);
bestAdapVecIndices = zeros(numberOfSubblocks,1);
pastLpcExcite = zeros(maxDelay,1);
bestAdapExcite = zeros(blockLength,1);
bestStochExcite = zeros(blockLength,1);
invLpcAdapExcite = zeros(blockLength,1);
exciteBlock = zeros(blockLength,1);

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

% Vector 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 );


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % LPC
  lpcInputBlock = xBlock;

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

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

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Analysis-by-Synthesis: 

  %
  % search for best adaptive excitation 
  %

  % take old filter states into account 
  zeroOutput = filter(1, [1;-aQ], zeroInput, zInvLpc);
  targetAdap = xBlock - zeroOutput; 

  [bestAdapGainIndices, bestAdapVecIndices] = ... 
      clltp( targetAdap, aQ, subblockLength, minDelay, adapGainCB, pastLpcExcite );

  bestAdapExcite = cradapex( subblockLength, ...
      bestAdapGainIndices, bestAdapVecIndices, adapGainCB, pastLpcExcite );

  %
  % search for best stochastic excitation
  %

  % take old filter states and adaptive excitation into account
  invLpcAdapExcite = filter(1, [1;-aQ], bestAdapExcite);
  targetStoch = targetAdap - invLpcAdapExcite; 

  % calculate all possible synthezized speech vectors 
  for l=1:stochVecCBSize,
    synFilterOutSignals(:,l) = filter( 1, [1;-aQ], stochVecCB(:,l) );
  end

  % gain-shape vector quantization of target   
  [bestStochGainIndex, bestStochVecIndex] = gsvecq( targetStoch, stochGainCB, synFilterOutSignals );


  % calculate new lpc filter states for next block 
  bestStochExcite = stochGainCB(bestStochGainIndex) * stochVecCB(:,bestStochVecIndex);
  exciteBlock = bestAdapExcite + bestStochExcite;
  [dummy, zInvLpc] = filter(1, [1;-aQ], exciteBlock, zInvLpc);

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

  % calculate new pastExcitation for next block
  pastLpcExciteExt = [pastLpcExcite; exciteBlock];
  pastLpcExcite = pastLpcExciteExt( blockLength+1 : maxDelay+blockLength ); 


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % simulate transmission of parameters via channel to decoder
  aChannel(:,i) = aQ;
  adapGainChannel(:,i) = bestAdapGainIndices;
  adapVecChannel(:,i) = bestAdapVecIndices;
  stochGainChannel(i) = bestStochGainIndex; 
  stochVecChannel(i) = bestStochVecIndex;

  blockIndex = blockIndex + blockLength;
end  % coder main loop


%
% THIS IS DECODER 
%

fprintf(1,'\nDecoding')

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

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

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


% decoder main loop
blockIndex = 1 : blockLength;

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

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % create excitation  

  bestAdapGainIndices = adapGainChannel(:,i);
  bestAdapVecIndices = adapVecChannel(:,i);
  bestAdapExcite = cradapex( subblockLength, ...
      bestAdapGainIndices, bestAdapVecIndices, adapGainCB, pastLpcExcite );

  bestStochGainIndex = stochGainChannel(i);
  bestStochVecIndex = stochVecChannel(i);
  bestStochExcite = stochGainCB(bestStochGainIndex) * stochVecCB(:,bestStochVecIndex);

  exciteBlock = bestAdapExcite + bestStochExcite;


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Inverse LPC (vocal tract filter)
  a = aChannel(:,i);
  [invLpcOutputBlock, zInvLpc] = filter(1, [1;-a], exciteBlock, zInvLpc);


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


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % calculate new pastExcitation for next looping
  pastLpcExciteExt = [pastLpcExcite; exciteBlock];
  pastLpcExcite = pastLpcExciteExt( blockLength+1 : maxDelay+blockLength ); 

  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')