cal: add functions for capsule calibration, fix gain normalization in case of negative values of power for S/T, adapt wobble fix to new version of findpeaks

calibration/src/sphear_calibration_common.m

@@ -81,6 +81,35 @@ function [AF, FS] = sphear_af_read(measurements, raw = false);
   end
 end
 
+function [AF, FS] = sphear_af_read_capsule(measurements, raw = false);
+  global DATA;
+  global CURRENT;
+  if (raw)
+    type = "";
+  else
+    type = "cal-";
+  end
+  if (exist([DATA CURRENT "af/capsule_u-cal-001-1.wav"]))
+    %% assume it is an eight capsule microphone
+    capnames = cellstr(['1'; '2'; '3'; '4'; '5'; '6'; '7'; '8']);
+    capsules = 1;
+  else
+    %% otherwise a tetrahedral microphone
+    capnames = cellstr(['LF'; 'RF'; 'LB'; 'RB']);
+    capsules = 4;
+  end
+  %% read all measurements
+  for m = measurements
+    %% read all capsules
+    for c = 1:capsules
+      %% read in calibrated capsule measurement
+      printf("m:%d -> c:%d\n", m, c);
+      [A, FS] = wavread([DATA CURRENT "af/capsule_u-" type sprintf("%03d", m) "-" char(capnames(c)) ".wav"]);
+      AF(m, c, :) = A;
+    end
+  end
+end
+
 %%
 %% Calculate power for signal in a frequency range
 %%
@@ -403,11 +432,55 @@ function [W X Y Z S T U V] = sphear_bf_extract_filter_shapes(B, measurements)
   end
 end
 
+%%
+%% Extract average shapes of capsule filters
+%%
+%% (using for now only the horizontal plane measurements)
+
+function [SHAPE] = sphear_extract_capsule_filter_shapes(A)
+  %% for each capsule
+  for c = (1:size(A,2))
+    %% get indexes of measurements to average for each capsule
+    mlll = mod(16 - (c * 2) + 0, 16) + 1;
+    mll = mod(16 - (c * 2) + 1, 16) + 1;
+    ml = mod(16 - (c * 2) + 2, 16) + 1;
+    m = mod(16 - (c * 2) + 3, 16) + 1;
+    mr = mod(16 - (c * 2) + 4, 16) + 1;
+    mrr = mod(16 - (c * 2) + 5, 16) + 1;
+    mrrr = mod(16 - (c * 2) + 6, 16) + 1;
+    for n = (1:5)
+      mm(n,:) = mod(16 - (c * 2) - n + 2, 16) + 1;
+      mm(n + 5,:) = mod(16 - (c * 2) + n + 2, 16) + 1;
+    end
+    printf("c:%d -> %d %d %d %d %d %d %d\n", c, mlll, mll, ml, m, mr, mrr, mrrr);
+    SHAPE(c,:) = (1.0 ./ reshape(mean(A([mll,ml,m,mr,mrr],c,:)), size(A,3), 1));
+    %% average all directions
+    %% SHAPE(c,:) = (1.0 ./ reshape(mean(A(:,c,:)), size(A,3), 1));
+  end
+end
+
+%%
+%% Filter signals through capsule correction filters
+%%
+
+function [AFF] = sphear_filter_capsule_signals(AF, FIR, fft_size);
+  measurements = size(AF, 1);
+  capsules = size(AF, 2);
+  samples = size(AF,3);
+  AFF = zeros(measurements, capsules, samples);
+  for m = (1:measurements)
+    for c = (1:capsules)
+      %% we add trailing zeros to let the filter ring
+      AFF(m,c,:) = fftfilt(reshape(FIR(c,:), 1, [])', cat(1, reshape(AF(m,c,:), 1, [])', zeros(2*size(FIR(1,:))(2),1))', fft_size)(1:samples);
+    end
+  end
+end
+
 %%
 %% Calculate average power for each B format component and normalize A2B matrix
 %%
 
-function [A2B XY_gain, S_gain, T_gain, U_gain, V_gain] = sphear_bf_normalize_gain(A2B, B, FS, fmin, fmax, fft_size)
+function [A2B XY_gain, S_gain, T_gain, U_gain, V_gain, B_power] = sphear_bf_normalize_gain(A2B, B, FS, fmin, fmax, fft_size)
   do_Z = false;
   if (size(B,1)/3 > 8)
     %% we assume we have three sets of measurements and can do Z averaging
@@ -434,11 +507,10 @@ function [A2B XY_gain, S_gain, T_gain, U_gain, V_gain] = sphear_bf_normalize_gai
     U_gain = U_avg/W_avg;
     V_avg = mean(mean(B_power([3,7,11,15], 8)));
     V_gain = V_avg/W_avg;
-    if (mean(mean(A2B(:,[5,6]))) > 0)
+    if (abs(mean(mean(A2B(:,[5,6])))) > 0)
       %% calculate peak value of s/t at measured elevation
       el = 43.0;
       s_t_val = 1.0 * sin((2*el)/180*pi);
-      s_t_val
       %% ST, peaks of lobes at +-elevation
       S_avg = mean([mean(B_power([1,9] + 16, 5)), mean(B_power([1,9] + 32, 5))]);
       S_gain = (S_avg/W_avg) * s_t_val;
@@ -460,7 +532,8 @@ function [A2B XY_gain, S_gain, T_gain, U_gain, V_gain] = sphear_bf_normalize_gai
     A2B(:,[2,3,4]) = A2B(:,[2,3,4]) / XY_gain;
   end
   if (size(A2B, 2) > 4)
-    if (mean(mean(A2B(:,[5,6]))) > 0)
+    if (abs(mean(mean(A2B(:,[5,6])))) > 0)
+      printf("correct ST %f %f\n", S_gain, T_gain);
       A2B(:,[5]) = A2B(:,[5]) / S_gain;
       A2B(:,[6]) = A2B(:,[6]) / T_gain;
     end
@@ -723,12 +796,25 @@ end
 %%%
 
 function [FIXED] = sphear_fix_wobble(SIG, DELAY);
+  sinc_size = 255;
   min_delay = min(DELAY);
   DELAY = DELAY - min_delay;
   max_delay = max(DELAY);
   for m = [1:size(DELAY, 1)]
-    printf("m: %d, delay: %f\n", m, DELAY(m));
-    FIXED(m,:,:) = sphear_delay_signal(SIG(m,:,:), max_delay - DELAY(m), 63, 8192);
+    %% delay = max_delay - DELAY(m);
+    delay = DELAY(m);
+
+%%    if (delay >= 1)
+%%      i_delay = floor(delay);
+%%      f_delay = mod(delay, 1.0);
+%%      printf("m: %d, delay: %f, integer %f, frac %f\n", m, delay, i_delay, f_delay);
+%%      FIXED(m,:,:) = zeros(size(SIG, 3));
+%%      FIXED(m,:,[i_delay:end]) = SIG(m,:,[1:end - i_delay])
+%%      FIXED(m,:,:) = sphear_delay_signal(FIXED(m,:,:), f_delay, 255, 8192);
+%%    else
+      printf("m: %d, delay: %f\n", m, delay);
+      FIXED(m,:,:) = sphear_delay_signal(SIG(m,:,:), delay, 255, 8192);
+%%    end
   end
 end