% Program to plot a signal's FFT spectrum, and the spectrum produced by
% a linear prediction model.
%

% file : /home/sb/sb1/speech_assign/assign5/lpc_hw.m
%

% main program starts here
%
sample_freq = 8000;
norm        = 3.0518e-05;
pi          = 3.1415927;
log_2_factor= 3.3219281;
msec_to_sec = 0.001;
bytes_per_sample = 2;

% get the parameters from the command line
%
frame_size = 20  % msec;
win_size   = 30  % msec;
win_center_time = 15 % msec;
num_chans  = 1;
preemp_value = 0.95 %input('   Preemphasis value  ');
lp_order    = 16   %input('   lp_order           ');
dft_order   = 2400 % (win_size * sample_freq * msec_to_sec);
win_type = 0  % (CLPC_RCTGL_WIN) 1 (CLPC_HAMMING_WIN)it should be a input value;

% calculate the number of samples per frame and per window
%
samples_per_win = (win_size * msec_to_sec * sample_freq);
win_center_sample = (win_center_time * msec_to_sec * sample_freq);

% open the audio file
%
[fid, message] = fopen('test.dat','r');

% set file pointers to required location in input file
%
[B, COUNT]= fread(fid,inf,'int16');

% compute the length of file in bytes
%
status = fseek(fid, 0, 1);
total_bytes = ftell(fid);

% end file sample 
%
end_file_sample = (total_bytes / (num_chans * bytes_per_sample));

% start and end sample of a window
%
offset_samples = round(samples_per_win/2);
start_sample   = round(win_center_sample - offset_samples);
end_sample     = round(win_center_sample + offset_samples);

% put this process in a loop
%
while (end_sample < end_file_sample) 

% declare and initialize buffer to hold 
% window of signal
%
num_samples = num_chans * samples_per_win;

% read in the signal center at win_center_time
%
new_start = start_sample;
new_end = end_sample;
num_samples = 1;

if (start_sample < 0)
for i=start_sample:1:1,
signal_window(num_samples) =0;
num_samples = num_samples + 1;
end
new_start =0;
end

if (end_sample > end_file_sample)
new_end = end_file_sample;
end

if(new_start > new_end)
fprintf(1,'start sample is greater than end sample\n');
end

% open the audio file
%
fid = fopen('test.dat','r');

% seek to the point in the file we want to read the
% data
fseek(fid, new_start * bytes_per_sample * num_chans, 0);
samples_read = (new_end - new_start) * num_chans;

% read these values 
%
[C, count] = fread(fid, samples_read, 'int16');

% convert the short int to floating point values
%
for i=1:samples_read,
signal_window(num_samples) = C(i);
num_samples = num_samples + 1;
end;

% pad it with 0s if window ends after file ends
%
%if (end_sample > end_file_sample) 
%  for i=end_file_sample:end_sample,
%	 signal_window(num_samples) = 0;
%num_samples = num_samples + 1;
%  end
%end

% debias the signal
% 
num_samples = samples_per_win * num_chans;

% compute the average signal value and subtract 
% from each original value
%
signal_avg = 0;

for i=1:num_chans,
for j=i:num_chans:num_samples,
signal_avg = signal_avg + signal_window(j);
end
signal_avg = signal_avg / samples_per_win;

% subtract this from sample values
%
for k=i:num_chans:num_samples,
signal_window(k) = signal_window(k) - signal_avg;
end
end

% iniitialize the pre_sig_val
%
for i=1:num_chans,
prev_sig_val(i) = 0;
end

% pre-emphasize the signal
%
%function pre_emphasize(prev_sig_val, samples_per_win, num_chans,preemp_val);
num_samples = samples_per_win * num_chans;

for i=1:num_samples,
temp_window(i) = signal_window(i);
end

% calculate the preemphasis signal
%
for i=1:num_chans,
signal_window(i) = temp_window(i) - preemp_value * prev_sig_val(i);
end

for i=1:num_chans,
for j=2:num_samples,
signal_window(j) = temp_window(j) - preemp_value * temp_window(j-num_chans);
end
end

for i=1:num_chans,
prev_sig_val(i) = signal_window(num_samples - num_chans + i);
end

% window the signal
%
%function window_signal(samples_per_win, num_chans, win_type);
num_samples = samples_per_win * num_chans;

if (win_type == 0)                  
for i=1:num_chans,
for j=i:num_chans:num_samples,
signal_window(j) = signal_window(j) * 1;
end
end


elseif (win_type == 1)
for i=1:num_chans,
for j=i:num_chans:num_samples,
signal_window(sample) = signal_window(sample) * (0.54 - 0.46 * cos(2 * pi * (i/num_chans)/(samples_per_win -1)));
end
end

else
fprintf(1,'Improper window type \n');
end

% offset for next window to be read
%
start_sample = start_sample + offset_samples;
end_sample = end_sample + offset_samples;

% while loop ends here
%
end

fprintf(1,'computing dft of spectrum \n');
% compute the DFT of the spectrum in dB
%
dft_samples = dft_order * num_chans;
signal_sample = samples_per_win * num_chans;

% pad signal with 0s if dft order is larger than
% window
for i=1:dft_samples,
if (i < signal_sample)
temp_window(i) = signal_window(i);
else
temp_window(i) =0;
end
end

for i=1:num_chans,
for j=1:dft_order,
twiddle_factor_rl = 0;
twiddle_factor_img =0;
samp =0;
for k=i:num_chans:dft_samples,
expn = 2.0 * pi * samp * j/dft_order;
twiddle_factor_rl = twiddle_factor_rl + temp_window(k) * cos(expn);
twiddle_factor_img = twiddle_factor_img - temp_window(k) * sin(expn);
samp = samp +1;
end
dft_spect_coeffs(i,j) = sqrt((twiddle_factor_rl * twiddle_factor_rl) + (twiddle_factor_img * twiddle_factor_img));
end
end

% open the file for writing the output
%
fid = fopen('output.text','w');

% print out the dft of signal to 
% stdout
%
num_samples = dft_order * num_chans;
freq_scaling = sample_freq / dft_order;

for i=1:num_chans,
for j=1:num_samples,
fprintf(fid, '%.4f  %.4f \n', freq_scaling * j, 20*log10(dft_spect_coeffs(i,j)));
end
end
fprintf(fid, '\n');

% compute the LPC spectrum in dB
%
%lp_signal = 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
lp_error = autocorr_coeffs(i,1);

% find the predictor coefficients using Levinson-Durbin recursion
%
refl_coeffs(j) = 0;
for j=2:lp_order,

% 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(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(i,j) = lp_signal(i,j) + sqrt(lp_gain(i) / (factor_rl * factor_rl + factor_img * factor_img));
end
end

% print out the lpc of signal to 
% stdout
%
num_samples = samples_per_win * num_chans;
freq_scaling = sample_freq / samples_per_win;

for i=1:num_chans,
for j=1:num_samples,
fprintf(fid, '%.4f %.4f \n', freq_scaling * j, 20*log10(lp_signal(i,j)));
end
end

%end
fclose(fid);

%
% file ends here































