% This function computes the lp of a spectrum
%
function lp_signal_fun = calculate_lp(signal_window, samples_per_win, num_chans, lp_order)

% compute the autocorrelation coefficients
%
autocorr_coeffs = (lp_order + 1);

% total number of samples for all channels
%
num_samples = (samples_per_win * num_chans);

for i=1:num_chans,
for j=1:autocorr_coeffs,
autocorr_coeffs(i,j)= 0;
sum =0;
for k=i:num_chans:num_samples - j * num_chans,
autocorr_coeffs(i,j) = autocorr_coeffs(i,j) + signal_window(k) * signal_window(k + j * num_chans);
sum = sum + autocorr_coeffs(i,j);
end 
end
end

% compute the lp predictor coefficients
%
num_coeffs = (lp_order + 1);
prev_pred_coeffs(num_coeffs) = 0;
for i=1:num_chans,

%initalize the error term


% find the predictor coefficients using Levinson-Durbin recursion
%
refl_coeffs(j) = 0;
for j=2:lp_order,
lp_error = sum;
% calculate and update the reflection coeffts
%
for k=1:j-1
  refl_coeffs(j) = refl_coeffs(j) - prev_pred_coeffs(k) * autocorr_coeffs(i,(j-k));
end

refl_coeffs(j) = refl_coeffs(j) - autocorr_coeffs(i,j);
refl_coeffs(j) = refl_coeffs(j) / lp_error;

% calculate and update the predictor coefficients
%
pred_coeffs(i,j) = refl_coeffs(j);
for k=1:j-1,
  pred_coeffs(i,k) = prev_pred_coeffs(k) + refl_coeffs(j) * prev_pred_coeffs(j-k);
end

% calculate and update the error
%
lp_error = lp_error * 1.0 - refl_coeffs(j) * refl_coeffs(j);

% update the temporary predictor coeffs
%
for k=1:j,
  prev_pred_coeffs(k) = pred_coeffs(i,k);
end
end

% invert the prediction coefficients
%
for j=1:lp_order,
  pred_coeffs(i,j) = pred_coeffs(i,j) * -1;
end

% set the gain value
%
lp_gain(i) = lp_error;
end

% computes the spectrum of signal using the coefficients
%
for i=1:num_chans,
  for j=1:samples_per_win,
  lp_signal_fun(i,j) =0;
end
end

for i=1:num_chans,
  for j=1:samples_per_win,
  factor_rl =1;
factor_img =0;

for k=1:lp_order,
  expn = (( 2 * pi / samples_per_win) * mod((k * j),samples_per_win));
factor_rl = factor_rl - pred_coeffs(i,k) * cos(expn);


factor_img = factor_img + pred_coeffs(i,k) * sin(expn);
end
lp_signal_fun(i,j) = lp_signal_fun(i,j) + sqrt(lp_gain(i) / (factor_rl * factor_rl + factor_img * factor_img));
end
end
return ;









