%%
%% *SpHEAR Project
%%
%% sphear_capsule_power: return sorted array of power for capsule measurements
%%
%% usage:
%%   [R Fs] = sphear_af_power(measurements, fmin, fmax, N, Npre)
%%
%% measurements: list of measurements to consider
%% fmin, fmax: range of frequencies to consider
%% Npre: number of samples before the main impulse
%% N: total number of samples to analyze
%%
%% The name of capsule measurement files iscurrently hardwired into this
%% file, this needs to be fixed and turned into input parameters
%%
%% Copyright Fernando Lopez-Lezcano 2016
%% nando@ccrma.stanford.edu
%%
%%  This program is free software: you can redistribute it and/or modify
%%  it under the terms of the GNU General Public License as published by
%%  the Free Software Foundation, either version 3 of the License, or
%%  (at your option) any later version.
%%
%%  This program is distributed in the hope that it will be useful,
%%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%%  GNU General Public License for more details.
%%
%%  You should have received a copy of the GNU General Public License
%%  along with this program.  If not, see <http://www.gnu.org/licenses/>.

function [RS Fs] = sphear_capsule_power(measurements, fmin, fmax, N, Npre)
  pkg load signal;

  onames = cellstr(['00'; '90']);

  %% all submitted measurements
  i = 1;
  for m = measurements
    %% orientation
    for o = [1, 2]

      printf("%s\n", ["../cal/2016.07.26/export/capsule-" sprintf("%02.0f", m) "-" char(onames(o)) "-1.wav"]);
      
      [capsule, Fs] = wavread(["../cal/2016.07.26/export/capsule-" sprintf("%02.0f", m) "-"  char(onames(o)) "-1.wav"]);

      %% indexes in array for frequency bins
      imin = round((N/2) * fmin / (Fs/2));
      imax = round((N/2) * fmax / (Fs/2));

      %% find peakr for this measurement
      [peaks, locs] = findpeaks(abs(capsule));
      [~, index] = max(peaks);
      start = locs(index) - Npre;
      range = start:(start + N - 1);
      
      fdom = fft(capsule(range), N);
      
      %% return total power in dB of signal in frequency range
      R(i, o) = 20*log10(sqrt(sum(fdom(imin:imax).*conj(fdom(imin:imax)))/size(fdom(imin:imax),1)));
    end
    %% add 0 to 90 degree difference in dB (front to side ratio)
    R(i, 3) = R(i, 1) - R(i, 2);
    if (i > 1)
      R(i, 4) = R(i, 3) - R(i - 1, 3);
    else
      R(i, 4) = 0;
    end
    R(i, 5) = i;
    i = i + 1;
  end
  avg = mean(R(:,1));
  R(:,1) = R(:,1) - avg;
  R(:,2) = R(:,2) - avg;
  %% sort the results by directivity
  [RS, I] = sort(R'(3,:));
  RS = R'(:,I)';
  %% calculate the derivative of directivity
  i = 1;
  for m = measurements
    if (i > 1)
      RS(i, 4) = RS(i, 3) - RS(i - 1, 3);
    else
      RS(i, 4) = 0;
    end
    %% add an index for plots
    RS(i, 6) = i;
    i = i + 1;
  end
end