function y = celp( x )
% function y = celp( x )
%
% Simulation of pure CELP Codec 
%
% 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:    03.01.2002 
% Contact: www.markus-hauenstein.de
%
% This program is copyrighted.

test = 0;

fprintf(1,'\nSimulation of pure CELP')
fprintf(1,'\nInitializing')

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

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

% size of stochastic gain codebook
stochGainCBSize = 64;

% size of stochastic Vector codebook
stochVecCBSize = 128;


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

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

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

% create normalized stochastic Vector codebook
stochCBAmp = 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 = stochCBAmp * stochVecCB;


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


% format Input signal
x=x(:); % make column Vec
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);
stochVecChannel = zeros(1, numberOfBlocks);
stochGainChannel = zeros(1, numberOfBlocks);


%
% THIS IS CODER 
%

fprintf(1,'\nCoding ')

% initialize coder variables

% lpc 
zInvLpc = zeros(lpcOrder,1);

% Vector quantization of LPC/LTP residual - helper variables for Vector quantization
synFilterOutSignals = zeros( blockLength, stochVecCBSize );
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
  lpcInputBlock = xBlock;

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

  % simulate quantization of LPC parameters   
  %aQ = larpredq( a, larMaxQ, predQuantBits ); 
  aQ = lsfpredq( a, predQuantBits ); 


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

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

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

  % Vector quantization of target   
  [bestStochGainIndex, bestStochVecIndex] = gsvecq( target, stochGainCB, synFilterOutSignals );


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % calculate lpc filter states for next looping 
  invLpcInputBlock = stochGainCB(bestStochGainIndex) * stochVecCB(:,bestStochVecIndex);
  [dummy, zInvLpc] = filter(1, [1;-aQ], invLpcInputBlock, zInvLpc);


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % simulate transmission of parameters via channel to decoder
  aChannel(:,i) = aQ;
  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 lpc
zInvLpc = zeros(lpcOrder,1);

% decoder main loop
blockIndex = 1 : blockLength;

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

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Inverse LPC
  invLpcInputBlock = stochGainCB(stochGainChannel(i)) * stochVecCB(:,stochVecChannel(i));
  a = aChannel(:,i);
  [invLpcOutputBlock, zInvLpc] = filter(1, [1;-a], invLpcInputBlock, zInvLpc);


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % construct output signal
  yBlock = invLpcOutputBlock; 
  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')