% Program to compute the signal to noise ratio of a raw
% speech file
%

% file : /home/sb/sb1/speech_assign/assign4/snr_calculate.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 = input('best frame size  (in msec)      ');  % 20 msec;
win_size   = input('best window size (in msec)      '); % 30 msec;
num_chans  = 1  % input values ;
noise_threshold = input('noise threshold %') ;
signal_threshold = input('signal threshold %');

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

% open the audio file
%
input_file = input('Enter the data file \n','s');
[fid, message] = fopen(input_file,'r');

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

% calculate the end file sample
%

% rewind the audio file
%
status = fseek(fid, 0, -1);

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

% calculate the number of frames in a file
%
total_frames = (round(total_bytes/(samples_per_frame * num_chans * bytes_per_sample)) + 1);

if (total_frames == 0)
     fprintf(1,' total frames is zero \n');
end

num_energy =1;

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

% counter for frame 
%
frame = 1;

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

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

% compute window of a signal
%
num_samples = num_chans * samples_per_win;

% read the signal frame by frame 
%
new_start = start_sample;
new_end   = end_sample;
num_samples = 1;

if (start_sample < 0 )
% pad with 0s if window start before file starts
%
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

fid = fopen('test4.dat','r');

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

[B,count] = fread(fid, samples_read, 'int16');

if (count ~= samples_read)
end

for i=1:samples_read,
signal_window(num_samples) = B(i);
num_samples = num_samples + 1;
end

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

% perform all the operations to process the signal data     
% declare the buffer and initialize all the items
%
for i=1:num_chans,
signal_prev_frame(i) = 0;
end

% perform pre_emphasizing the signal
%
num_samples = samples_per_win * num_chans;
for i=1:num_chans,
pre_emphasized_signal(i) = signal_window(i) - 0.95 * signal_prev_frame(i);
end

for j=1:num_chans,
  for i=num_chans + j:num_chans:num_samples,
  pre_emphasized_signal(i) = signal_window(i) - 0.95 * signal_window(i-num_chans);
end
end

for j=1:num_chans,
  signal_prev_frame(j) = signal_window(num_samples - num_chans + j);
end

% perform hamming window operation
%
for j=1:num_chans,
  for i=1:num_chans:num_samples,
  hamming_window_signal(i) = signal_window(i) * (0.54 -0.46 * cos(2*3.14 * (i/num_chans) / (samples_per_win - 1)));
end
end

% perform the computation of energy
%
for j=1:num_chans,
  energy_signal(j) =0;
for i=j:num_samples,
  energy_signal(j) = energy_signal(j) + signal_window(i) * signal_window(i);
end
energy_signal(j) = energy_signal(j) / samples_per_win;
end

% perform add energy 
%
for j=1:num_chans,
energy_buffer(j,num_energy) = energy_signal(j);
end
num_energy = num_energy + 1;

% offset for next frame to be read
%
start_sample = start_sample + samples_per_frame;
end_sample = end_sample + samples_per_frame;
frame = frame +1 ;

% while loop ends here
%
end

% perform the pdf and cdf computation
%

% compute the histogram of the energy of signal of each channel
%
max = 100000;
min = -100000;

for i=1:num_energy-1,
  if (energy_buffer(i) > max)
     max_value(i) = energy_buffer(i);
     else
max_value(i) = max;
end

if (energy_buffer(i) < min)
     min_value(i) = energy_buffer(i)
     else
min_value(i) = min;
end
end   

for i=1:num_chans,
  bin_width(i) = (max_value(i) - min_value(i)) / (10000);
end

% initialize the bins
     %
     for j=1:num_chans,
       for i=1:10000,
       bin_freqs(j,i) = 0;
bin_avgs(j,i) = 0;
end
end

% determine what goes to which bin
%
for j=1:num_chans,
  for i=1:num_energy-1,
  
  % limit the data
%
num = energy_buffer(j,i);
if (num < min_value(j))
    num = min_value(j);
    
    elseif (num > max_value(j))
    num = max_value(j);
end

% calculate which bin the data should go and increment the count
%
index = round((num - min_value(j)) / bin_width(j));
bin_freqs(j,index) = bin_freqs(j,index) + 1;
bin_avgs(j,index) = bin_avgs(j,index) + num;
end
end

for j=1:num_chans,
for i=1:10000,

if (bin_freqs(j,i) == 0)
bin_avgs(j,i) = (bin_width(j) * i - bin_width(j)) / 2;
else
bin_avgs(j,i) = bin_avgs(j,i) / bin_freqs(j,i);
end
end
end
% compute the pdf of the energy of signal
%
for j=1:num_chans,
for i=1: 10000,
pdf(j,i) = bin_freqs(j,i) / num_energy;
end
end

% compute the cdf of the energy of signal
%
for j=1:num_chans,
cdf(j,1) = pdf(j,1);

% loop over all bins
%
for i=2:10000,
cdf(j,i) = pdf(j,i) + cdf(j,i-1);
end
end

% compute the snr value for each channel
%

% set the flag values
%
noise_flag = 0 % false;
signal_flag = 0 % false;

for j=1:num_chans,
for i=1:10000,
if (cdf(j,i)>(noise_threshold)) & (noise_flag==0)
noise_level = bin_avgs(j,i);
noise_flag = 1;

elseif ((cdf(j,i) > signal_threshold) & (signal_flag==0))
signal_level = bin_avgs(j,i);
signal_flag = 1;
end
end
snr_values(j) = 10 * log10((signal_level - noise_level) / noise_level);
noise_flag =0;
noise_flag =0;
end

% output the snr values for each channel
%
for i=1: num_chans,
fprintf(1,'snr value for this audio file is : %.6f dB \n', snr_values(i));
end

% close the file
%
fclose(fid);






