function y = lpcvocod( x )
% function y = lpcvocod( x )
%
% Simulation of simple LPC Vocoder 
%
% 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:    18.05.2002 
% Contact: www.markus-hauenstein.de
%
% This program is copyrighted.

test = 0;

fprintf(1,'\nSimulation of simple LPC Vocoder');
fprintf(1,'\nInitializing');


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

% general parameters
blockLength = 160;

% gain quantization
ampQuantBits = 5;

% ltp (number of samples)
minDelay = 20;
maxDelay = 147;
akfBlockLength = 40;

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


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

% gain quantization
maxAmp = 32000;

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


% calculate and display required bitrate
sampleFrequency = 8000;
stateBitRate = 2 * sampleFrequency / blockLength; 
ampBitrate = ampQuantBits * sampleFrequency / blockLength;
lpcBitRate = sum( predQuantBits ) * sampleFrequency / blockLength;
ltpDelayBitRate = ceil( log(maxDelay-minDelay)/log(2) ) * sampleFrequency / blockLength;
bitRate = stateBitRate + ampBitrate + lpcBitRate + ltpDelayBitRate;
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 
stateChannel = zeros(1, numberOfBlocks);
ampChannel = zeros(1, numberOfBlocks);
aChannel = zeros( lpcOrder, numberOfBlocks);
ltpDelayChannel = zeros(1, numberOfBlocks);

% initialize classifier parameters
globalMaxAmp = 1024;
silenceMaxLevel = -60;
upperVoicedLimit = 0.3;

% initialize variables for encoding frame types 
silenceFrame = 1;
voicedFrame = 2;
voicelessFrame = 3;

% initialize variables for encoding coder states
silence = 1;
voiced = 2;
voiceless = 3;
wasVoiced = 4;
wasVoiceless = 5;


%
% THIS IS CODER 
%

fprintf(1,'\nCoding'); 
state = silence;

% coder main loop
blockIndex = 1 : blockLength;

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

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

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % silence / voiced / voiceless decision

  % classification of current block

  globalMaxAmp = max( globalMaxAmp, max(abs(xBlock)) );
  xAmp = sqrt( xBlock' * xBlock ) / blockLength;
  level = 20*log10( xAmp / globalMaxAmp );

  if level < silenceMaxLevel
    current = silenceFrame;
  else
    % check if voiced or voiceless frame
    signum = sign( xBlock );
    relZeroCrossRate = ...
        sum( 0.5 * abs( signum(1:blockLength-1) - signum(2:blockLength) ) ) / (blockLength-1);
    if relZeroCrossRate <= upperVoicedLimit 
      current = voicedFrame;
    else
      current = voicelessFrame;
    end
  end

  % state machine which takes history into account  
  % 
  if     state == silence      & current == silenceFrame    state = silence;
  elseif state == silence      & current == voicedFrame     state = voiced;
  elseif state == silence      & current == voicelessFrame  state = voiceless;
  elseif state == voiced       & current == silenceFrame    state = wasVoiced;
  elseif state == voiced       & current == voicedFrame     state = voiced;
  elseif state == voiced       & current == voicelessFrame  state = voiceless;
  elseif state == voiceless    & current == silenceFrame    state = wasVoiceless;
  elseif state == voiceless    & current == voicedFrame     state = voiced;
  elseif state == voiceless    & current == voicelessFrame  state = voiceless;
  elseif state == wasVoiced    & current == silenceFrame    state = silence;
  elseif state == wasVoiced    & current == voicedFrame     state = voiced;
  elseif state == wasVoiced    & current == voicelessFrame  state = voiceless;
  elseif state == wasVoiceless & current == silenceFrame    state = silence;
  elseif state == wasVoiceless & current == voicedFrame     state = voiced;
  elseif state == wasVoiceless & current == voicelessFrame  state = voiceless;
  else disp('error in state machine');
  end       

  %state

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%%%%%%%%%%%%  Determine now paramenters, depending on state  %%%%%%%%%%%%%%%%
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %
  if state == silence
      stateChannel(:,i) = silence;
  end


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %
  if state == voiceless | state == wasVoiceless 

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % calculate and quantize amplitude of input signal
      %
      xAmp = sqrt( xBlock' * xBlock / blockLength );
      xAmpQ = muquant( xAmp, maxAmp, ampQuantBits );


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


      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % simulate transmission of parameters via channel to decoder

      stateChannel(:,i) = voiceless;
      ampChannel(:,i) = xAmpQ;
      aChannel(:,i) = aQ;
  end       


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %
  if state == voiced | state == wasVoiced

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % calculate and quantize amplitude of input signal
      %
      xAmp = sqrt( xBlock' * xBlock / blockLength );
      xAmpQ = muquant( xAmp, maxAmp, ampQuantBits );


      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % lpc coefficients (inverse vocal tract filter)

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


      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % ltp lag - simple block autocorrelation approach

      blockAkf = conv( xBlock(akfBlockLength:-1:1), xBlock );
      blockAkf = blockAkf( akfBlockLength : length(blockAkf) - akfBlockLength + 1 );     

      searchArea = blockAkf( minDelay + 1 : min(maxDelay, length(blockAkf)) );
      [maximum, maximumIndex] = max( searchArea );
      ltpDelay = maximumIndex(1) + minDelay - 1;


      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % simulate transmission of parameters via channel to decoder

      stateChannel(:,i) = voiced;
      ampChannel(:,i) = xAmpQ;
      aChannel(:,i) = aQ;
      ltpDelayChannel(:,i) = ltpDelay;
  end

end  % coder main loop



%
% THIS IS DECODER 
%

fprintf(1,'\nDecoding')

% initialize decoder variables
y = zeros(numberOfBlocks*blockLength,1);
randn('seed',0);
inputNoise = randn( blockLength, 1 );
zInvLpc = zeros(1, lpcOrder );
overhang = 0;


% decoder main loop
blockIndex = 1 : blockLength;

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

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % construct excitation signal for vocal tract filter
  
  state = stateChannel(i);

  %%%% if silence frame
  %
  if state == silence
      excitation = zeros( blockLength, 1 );
      overhang = 0;
  end

  %%%% if voiceless frame
  %
  if state == voiceless
      excitation = inputNoise;
      overhang = 0;
  end

  %%%% if voiced frame
  %
  if state == voiced
      excitation = zeros( blockLength, 1 );
      ltpDelay = ltpDelayChannel(i);
      
      if overhang == 0
          position = 1;
      else
          position = max(1, ltpDelay - overhang); 
          % overhang considers last pulse position in last frame
      end     

      while position <= blockLength,
          excitation( position ) = 1; 
          overhang = blockLength - position;
          position = position + ltpDelay;
      end     
  end  

  %plot( excitation )
  %pause


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


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % scaling
  invLpcPower = invLpcOutputBlock' * invLpcOutputBlock / blockLength;
  if invLpcPower > 0
    gain = sqrt( ampChannel(i) * ampChannel(i) / invLpcPower );
  else
    gain = 0.0; 
  end
  
  yBlock = gain * invLpcOutputBlock;


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


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % prepare next looping
  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');