% file: compute_snr.m
%
% this function computes the SNR (signal to noise ratio) of a raw data file
%
% arguments:
%  window_duration	duration of the window in msec
%  frame_duration	duration of the frame in msec
%  sample_rate		rate of the input data in kHz
%  in_file		binary input file
%  emphasis_val		parameter used in pre-emphasis
%  percentile_a		lower percentile number for snr (<1)
%  percentile_b		higher percentile number for snr (<1)
%  bins			number of bins to use for the histogram
%
% return value: 
%  snr			signal to noise ratio
%
% output: 
%  plot			energy histogram
%  plot			energy cdf
%

function snr = compute_snr(window_duration, frame_duration, sample_rate, in_file, emphasis_val, percentile_a, percentile_b, bins)

% get the durations in msecs
%
window_duration = window_duration/1000;
frame_duration = frame_duration/1000;

% open the raw data file
%
fp = fopen(in_file,'r');

% determine the number of samples per window and per frame
%
num_samp_window = window_duration*sample_rate + 1;
num_samp_frame = frame_duration*sample_rate + 1;

% create Hamming Window to window the raw data
%
hamming_window = hamming(num_samp_window);

% read in the first frame of raw data (16-bit int),
% zeropad the front part of the window
%
[first_frame,count] = fread(fp,num_samp_frame,'int16');
if ( num_samp_window ~= num_samp_frame )
   zero_pad = [zeros(floor((num_samp_window-num_samp_frame)/2),1)];
   [first_win] = fread(fp,ceil((num_samp_window-num_samp_frame)/2),'int16');
   data = [zero_pad; first_frame; first_win];
else
   data = first_frame;
end

% initialize count for the energy
%
energy_index = 1;

while (count ~= 0)

   % pre-emphasize the data
   %
   pe_data = pre_emph(data,emphasis_val);

   % apply Hamming Window to the data
   %
   windowed_data = hamming_window .* transpose(pe_data);

   % compute the energy 
   %
   energy(energy_index) = compute_energy(windowed_data);

   % read in the next data to be processed
   %
   [new_data,count] = fread(fp,num_samp_frame-1,'int16');

   % shift the data by one frame
   %
   if (count ~= 0)
      data = cat(1,data(num_samp_frame:num_samp_window),new_data);

      % zero-pad if necessary for the last frame of data
      %
      if ( length(data) < length(hamming_window) )
         data(length(hamming_window)) = 0;
      end

      energy_index = energy_index + 1; 
   end

end

% close the raw data file
%
fclose(fp);

% compute the pdf of the energy
%
energy_pdf = hist(energy,bins);

% compute the cdf of the energy
%
% 
energy_cdf = cumsum(energy_pdf);

% plot the energy histogram
%
% energy_hist_plot = figure;
% plot(1:bins,energy_pdf);
% grid on;
% title('Energy Histogram');

% plot the energy cdf
%
% energy_cdf_plot = figure;
% plot(1:bins,energy_cdf);
% grid on;
% title('Energy Cdf');

% compute the values of the nominal noise level 
% and the nominal signal+noise level
%
[nominal_noise, nominal_sig_noise] = compute_noise_levels(energy,energy_cdf,percentile_a,percentile_b,bins);


% compute the SNR
%
snr = 10*log10((nominal_sig_noise-nominal_noise)/nominal_noise);


