cal: many changes to the support function to sort capsule measurements for matching them

calibration/src/sphear_capsule_power.m

@@ -1,3 +1,4 @@
+%% -*- mode: octave -*-
 %%
 %% *SpHEAR Project
 %%
@@ -14,7 +15,7 @@
 %% 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
+%% Copyright Fernando Lopez-Lezcano 2016-2018
 %% nando@ccrma.stanford.edu
 %%
 %%  This program is free software: you can redistribute it and/or modify
@@ -30,11 +31,10 @@
 %%  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)
+function [RS Fs] = sphear_capsule_power(measurements, fmin, fmax, Nfft, Npre, Npost)
   pkg load signal;
 
   capname = "EM200";
-  
   onames = cellstr(['001'; '002']);
 
   %% all submitted measurements
@@ -44,23 +44,27 @@ function [RS Fs] = sphear_capsule_power(measurements, fmin, fmax, N, Npre)
     for o = [1, 2]
 
       printf("%s\n", ["../data/current/" capname "/export/capsule-" sprintf("%02.0f", m) "-" char(onames(o)) "-1.wav"]);
-      
-      [capsule, Fs] = wavread(["../data/current/" capname "/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));
+      if (exist(["../data/current/" capname "/export/capsule-" sprintf("%02.0f", m) "-"  char(onames(o)) "-1.wav"]))
+	[capsule, Fs] = wavread(["../data/current/" capname "/export/capsule-" sprintf("%02.0f", m) "-"  char(onames(o)) "-1.wav"]);
 
-      %% 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)))/(imax - imin)));
+	%% indexes in array for frequency bins
+	imin = round((Nfft/2) * fmin / (Fs/2));
+	imax = round((Nfft/2) * fmax / (Fs/2));
+
+	%% find peak for this measurement
+	[peaks, locs] = findpeaks(abs(capsule));
+	[~, index] = max(peaks);
+	start = locs(index) - Npre;
+	range = start:(start + Npost);
+	fdom = fft(capsule(range), Nfft);
+
+	%% return total power in dB of signal in frequency range
+	R(i, o) = 20*log10(sqrt(sum(fdom(imin:imax).*conj(fdom(imin:imax)))/(imax - imin)));
+      else
+	printf("missing capsule %d\n", m);
+	R(i, o) = -100;
+      end
     end
     %% add 0 to 90 degree difference in dB (front to side ratio)
     R(i, 3) = R(i, 1) - R(i, 2);
@@ -72,13 +76,14 @@ function [RS Fs] = sphear_capsule_power(measurements, fmin, fmax, N, Npre)
     R(i, 5) = i;
     i = i + 1;
   end
+  printf("calculating results...\n");
   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, I] = sort(R'(3,:));
   %% sort the results by 0 degree gain
-  [RS, I] = sort(R'(1,:));
+  %%[RS, I] = sort(R'(1,:));
   RS = R'(:,I)';
   %% calculate the derivative of directivity
   i = 1;