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