add function to process and sort capsule measurements

sphear_capsule_power.m

@@ -0,0 +1,92 @@
+%%
+%% *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
+