add preliminary support for Octathingy calibration and second order horizontal Ambisonics

calibration/src/sphear_calibration_common.m

@@ -65,19 +65,29 @@ end
 %%
 %% (with an option to replace Z by a forward facing cardioid if the signal is B format)
 
-function [R] = sphear_signal_power(SIG, FS, fmin, fmax, fft_size, cardioid = 0)
+function [R] = sphear_signal_power(SIG, FS, fmin, fmax, fft_size, cardioid = 0, azim = 45)
   %% indexes in array for frequency bins
   imin = round((fft_size/2) * fmin / (FS/2));
   imax = round((fft_size/2) * fmax / (FS/2));
   measurements = size(SIG)(1);
-  signals = size(SIG)(2);
+  signals = size(SIG,2);
   %% calculate power in frequency range
   for m = (1:measurements)
     for s = (1:signals)
-      if (s == 4 && cardioid)
-	spectrum = fft((SIG(m, 1, :) + 1.0 * SIG(m, 2, :)) / 2, fft_size);
+      if (s == 3 && cardioid > 1)
+	spectrum = fft((SIG(m, 1, :) +
+			(1.0 * SIG(m, 2, :) * cos(azim/180*pi)) +
+			(1.0 * SIG(m, 3, :) * sin(azim/180*pi)) +
+			(1.0 * SIG(m, 5, :) * cos((2 * azim)/180*pi)) +
+			(1.0 * SIG(m, 6, :) * sin((2 * azim)/180*pi))) / 2.5, fft_size);
+	%% spectrum = fft(SIG(m, s, :), fft_size);
       else
-	spectrum = fft(SIG(m, s, :), fft_size);
+	if (s == 4 && cardioid)
+	  %% spectrum = fft((SIG(m, 1, :) + 1.0 * SIG(m, 2, :)) / 2, fft_size);
+	  spectrum = fft((SIG(m, 1, :) + 1.0 * SIG(m, 2, :) * cos(azim/180*pi) + 1.0 * SIG(m, 3, :) * sin(azim/180*pi)) / 2, fft_size);
+	else
+	  spectrum = fft(SIG(m, s, :), fft_size);
+	end
       end
       %% return total power of signal in frequency range
       R(m,s,:) = sqrt(sum(spectrum(imin:imax).*conj(spectrum(imin:imax)))/(imax - imin));
@@ -97,6 +107,18 @@ function [A2B PV COND ERROR] = sphear_calculate_a2b(R)
   COND = cond(A2B);
 end
 
+function [A2B PV COND ERROR] = sphear_octa_calculate_a2b(R, uv_factor)
+  [A2B PV ERROR] = A2B_octa_matrix(R);
+  %% normalize matrix
+  %
+  %% ignore anything above first order as that will be difference
+  %% microphones that will be low level in the reference frequency
+  %% range and will skew the mean
+  A2B = A2B ./ mean(mean(abs(A2B(:,1:4))));
+  %% return condition number
+  COND = cond(A2B);
+end
+
 %% Render a calibrated capsule impulse response through an A2B matrix
 %%
 %% m: number of measurement to render
@@ -116,11 +138,13 @@ function [BF, BFF] = sphear_af_to_bf(AF, A2B, FIR, fft_size);
   if exist("FIR")
     %% filter the B format signal through the four FIR filters
     for m = (1:measurements)
-      for b = 1:4
+      for b = 1:size(FIR, 1)
 	%% we add trailing zeros to let the filter ring
 	BFF(m,b,:) = fftfilt(reshape(FIR(b,:), 1, [])', cat(1, reshape(BF(m,b,:), 1, [])', zeros(2*size(FIR(1,:))(2),1))', fft_size);
       end
     end
+  else
+    BFF = BF;
   end
 end
 
@@ -145,7 +169,7 @@ function [SP, FC, FE] = sphear_signal_power_all(SIG, FS, fmin, fmax, frefmin, fr
   %% (measurement,signal,power)
   ref = sphear_signal_power(SIG, FS, frefmin, frefmax, fft_size);
   norm = mean(ref(:,:));
-  
+
   for i = (1:size(FC))
     %% (measurement,signal,band,power)
     SP(:,:,i) = sphear_signal_power(SIG, FS, FE(i), FE(i + 1), fft_size) ./ norm;
@@ -163,12 +187,84 @@ TETRA_DIAGONALS = [3, 7, 11, 15];
 OCTA_CARDINALS = [1, 3, 5, 7, 9, 11, 13, 15];
 OCTA_DIAGONALS = [2, 4, 6, 8, 10, 12, 14, 16];
 
-function [W X Y Z] = sphear_bf_extract_filter_shapes(B, measurements)
+%% measurements:
+%% 0 -> use defaults
+%% 1 -> use all measurements
+%% 2 -> use cardinal measurements
+%% 3 -> use diagonal measurements
+
+function [W X Y Z U V] = sphear_bf_extract_filter_shapes(B, measurements)
+  if (measurements == 0)
+    if (size(B, 2) > 4)
+      W_m = [1, 3, 5, 7, 9, 11, 13, 15];
+      X_m = [4, 5, 6, 12, 13, 14];
+      Y_m = [16, 1, 2, 8, 9, 10];
+      Z_m = W_m;
+      %%U_m = [1, 5, 9, 13];
+      %%V_m = [3, 7, 11, 15];
+      U_m = [1:16];
+      V_m = [1:16];
+    else
+      W_m = [1, 5, 9, 13];
+      X_m = [1, 9];
+      Y_m = [5, 13];
+      Z_m = W_m;
+      U_m = [1, 5, 9, 13];
+      V_m = [3, 7, 11, 15];
+    end
+  else
+    if (measurements == 1)
+      W_m = [1:16];
+      X_m = [1:16];
+      Y_m = [1:16];
+      Z_m = [1:16];
+      U_m = [1:16];
+      V_m = [1:16];
+    else
+      if (measurements == 2)
+	W_m = [1, 5, 9, 13];
+	X_m = [1, 5, 9, 13];
+	Y_m = [1, 5, 9, 13];
+	Z_m = [1, 5, 9, 13];
+      else
+      end
+    end
+  end
   %% extract average response shape for W
-  W = 1.0 ./ reshape(mean(B(measurements,1,:)), size(B(:,:,:))(3), 1);
+  W = 1.0 ./ reshape(mean(B(W_m,1,:)), size(B(:,:,:))(3), 1);
   %% extract average response shape for XYZ
-  X = 1.0 ./ reshape(mean(B(measurements,2,:)), size(B(:,:,:))(3), 1);
-  Y = 1.0 ./ reshape(mean(B(measurements,3,:)), size(B(:,:,:))(3), 1);
-  Z = 1.0 ./ reshape(mean(B(measurements,4,:)), size(B(:,:,:))(3), 1);
+  X = 1.0 ./ reshape(mean(B(X_m,2,:)), size(B(:,:,:))(3), 1);
+  Y = 1.0 ./ reshape(mean(B(Y_m,3,:)), size(B(:,:,:))(3), 1);
+  Z = 1.0 ./ reshape(mean(B(Z_m,4,:)), size(B(:,:,:))(3), 1);
+  if (size(B, 2) > 4)
+    U = 1.0 ./ reshape(mean(B(U_m,5,:)), size(B(:,:,:))(3), 1);
+    V = 1.0 ./ reshape(mean(B(V_m,6,:)), size(B(:,:,:))(3), 1);
+  end
 end
 
+%%
+%% Calculate average power for each B format component and normalize A2B matrix
+%%
+
+function [A2B XY_gain, UV_gain] = sphear_bf_normalize_gain(A2B, B, FS, fmin, fmax, fft_size)
+  %% (measurement,signal,band,power)
+  B_power = sphear_signal_power(B, FS, fmin, fmax, fft_size);
+  %% W, all measurements
+  W_avg = mean(B_power(:,1));
+  %% XY, first order lobes
+  XY_avg = mean([mean(B_power([1,9],2)), mean(B_power([5, 13], 3))]);
+  if (size(A2B, 1) > 4)
+    %% UV, second order lobes
+    UV_avg = mean([mean(B_power([1,5,9,13],5)), mean(B_power([3,7,11,15], 6))]);
+  end
+
+  %% now calculate the gain coefficients and normalize the A2B matrix
+  XY_gain = XY_avg/W_avg;
+  A2B(:,[2,3,4]) = A2B(:,[2,3,4]) / XY_gain;
+  if (size(A2B, 1) > 4)
+    UV_gain = UV_avg/W_avg;
+    A2B(:,[5,6]) = A2B(:,[5,6]) / UV_gain;
+  else
+    UV_gain = 1.0;
+  end
+end