calibration: add preliminary support for ST second order components for Octathingy calibration

calibration/src/A2B_matrix.m

@@ -35,13 +35,26 @@ function [A2B pv max_err] = A2B_matrix(Y, X, v_z_factor)
   X1 = [ ones(size(Y,1), 1) X ];
   pv = X1 \ Y;
   max_err = max (Y - (X1 * pv));
+  %% fill in Z projections if needed
+  if(sum((X1(:,4))==0) > 0)
+    pv = fill_v_z( pv );
+  end
   if (size(Y, 2) > 4)
-    %% overdetermined
-    pv = fill_v_z( pv )';
-    ipv = pinv(pv);
-    A2B = ([1 0 0 0 0 0; 0 1 0 0 0 0; 0 0 1 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 1] * ipv)';
+    ipv = pinv(pv');
+    if (size(pv, 1) > 6)
+      %% W X Y Z S T U V
+      A2B = ([1 0 0 0 0 0 0 0; 0 1 0 0 0 0 0 0; 0 0 1 0 0 0 0 0; 0 0 0 1 0 0 0 0;
+	      0 0 0 0 1 0 0 0; 0 0 0 0 0 1 0 0; 0 0 0 0 0 0 1 0; 0 0 0 0 0 0 0 1] * ipv)';
+    else
+      %% W X Y Z U V
+      A2BI = ([1 0 0 0 0 0; 0 1 0 0 0 0; 0 0 1 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 1] * ipv)';
+      %% fill in S/T components with zeros
+      A2B(:,[1,2,3,4]) = A2BI(:,[1,2,3,4]);
+      A2B(:,[5,6]) = zeros(8, 2);
+      A2B(:,[7,8]) = A2BI(:,[5,6]);
+    end
   else
-    pv = fill_v_z( pv );
+    %% W X Y Z
     A2B = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1] / pv;
   end
 end

calibration/src/sphear_bf_plot.m

@@ -47,12 +47,16 @@ function [BF] = sphear_bf_plot(BF, fmin, fmax, rmax, cardioid = 0)
   if (cardioid)
     bname = cellstr(['W'; 'X'; 'Y'; 'Cardioid'; 'U'; 'V']);
   else
-    bname = cellstr(['W'; 'X'; 'Y'; 'Z'; 'U'; 'V']);
+    if (size(BF, 2) > 6)
+      bname = cellstr(['W'; 'X'; 'Y'; 'Z'; 'S'; 'T'; 'U'; 'V']);
+    else
+      bname = cellstr(['W'; 'X'; 'Y'; 'Z'; 'U'; 'V']);
+    end
   end
 
   f = figure;
   for b = 1:size(BF, 2)
-    subplot(2, 3, b);
+    subplot(2, 4, b);
     hold on;
     grid off;
     axis square;

calibration/src/sphear_bf_plot_frequency_response.m

@@ -73,9 +73,11 @@ function [] = sphear_bf_plot_frequency_response(BF, frefmin, frefmax, fmin, fmax
     semilogx(fr, bfft(m,3,1:end/2), 'b;Y;', "linewidth", 1.5);
     semilogx(fr, bfft(m,4,1:end/2), 'g;Z;', "linewidth", 1.5);
     if (size(BF, 2) > 4)
-      semilogx(fr, bfft(m,5,1:end/2), 'r--;U;', "linewidth", 1.5);
-      semilogx(fr, bfft(m,6,1:end/2), 'b--;V;', "linewidth", 1.5);
-      l = legend('W', 'X', 'Y', 'Z', 'U', 'V');
+      semilogx(fr, bfft(m,5,1:end/2), 'c:;S;', "linewidth", 1.5);
+      semilogx(fr, bfft(m,6,1:end/2), 'm:;T;', "linewidth", 1.5);
+      semilogx(fr, bfft(m,7,1:end/2), 'r--;U;', "linewidth", 1.5);
+      semilogx(fr, bfft(m,8,1:end/2), 'b--;V;', "linewidth", 1.5);
+      l = legend('W', 'X', 'Y', 'Z', 'S', 'T', 'U', 'V');
     else
       l = legend('W', 'X', 'Y', 'Z');
     end

calibration/src/sphear_calibration_common.m

@@ -58,6 +58,27 @@ function [AF, FS] = sphear_af_read(measurements, raw = false);
       AF(m, c, :) = A;
     end
   end
+  %% read any up/down measurements
+  if (exist([DATA CURRENT "af/a-format-cal-up-001-5.wav"]))
+    for m = measurements
+      %% read all capsules
+      for c = 1:capsules
+	%% read in calibrated capsule measurement
+	[A, FS] = wavread([DATA CURRENT "af/a-format-" type sprintf("up-%03d", m) "-" char(capnames(c)) ".wav"]);
+	AF(m + 16, c, :) = A;
+      end
+    end
+  end
+  if (exist([DATA CURRENT "af/a-format-cal-down-001-5.wav"]))
+    for m = measurements
+      %% read all capsules
+      for c = 1:capsules
+	%% read in calibrated capsule measurement
+	[A, FS] = wavread([DATA CURRENT "af/a-format-" type sprintf("down-%03d", m) "-" char(capnames(c)) ".wav"]);
+	AF(m + 16 + 16, c, :) = A;
+      end
+    end
+  end
 end
 
 %%
@@ -98,21 +119,52 @@ function [R] = sphear_signal_power(SIG, FS, fmin, fmax, fft_size, cardioid = 0,
   end
 end
 
+%%
+%% Calculate power for capsule in a frequency range
+%%
+
+function [R] = sphear_capsule_power(SIG, FS, fmin, fmax, fft_size)
+  %% 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);
+  %% calculate power in frequency range
+  for m = (1:measurements)
+    for s = (1:signals)
+      spectrum = fft(SIG(m, s, :), fft_size);
+      %% return total power of signal in frequency range
+      R(m,s,:) = sqrt(sum(spectrum(imin:imax).*conj(spectrum(imin:imax)))/(imax - imin));
+    end
+  end
+end
+
 %%
 %% Calculate A to B matrix from matrix of power measurements
 %%
 
-function [A2B PV COND ERROR] = sphear_calculate_a2b(R)
-  [A2B PV ERROR] = A2B_matrix(R);
+function [A2B PV COND ERROR] = sphear_calculate_a2b(R, M)
+  if (exist("M"))
+    [A2B PV ERROR] = A2B_matrix(R, M);
+  else
+    [A2B PV ERROR] = A2B_matrix(R);
+  end
   %% normalize matrix
   %%
   %% normalize first order components of the matrix first
   A2B = A2B ./ mean(mean(abs(A2B(:,1:4))));
   if (size(A2B, 2) > 4)
-    %% U/V are difference microphones, so at the reference
+    %% S/T/U/V are difference microphones, so at the reference
     %% frequency range their signal is much lower than the
     %% first order components, so normalize separately
-    A2B(:,5:6) = A2B(:,5:6) ./ mean(mean(abs(A2B(:,5:6))));
+    ST_avg = mean(mean(abs(A2B(:,5:6))));
+    if (ST_avg != 0)
+      A2B(:,5:6) = A2B(:,5:6) ./ ST_avg;
+    end
+    UV_avg = mean(mean(abs(A2B(:,7:8))));
+    if (UV_avg != 0)
+      A2B(:,7:8) = A2B(:,7:8) ./ UV_avg;
+    end
   end
   %% return condition number
   COND = cond(A2B);
@@ -135,7 +187,7 @@ function [BF, BFF] = sphear_af_to_bf(AF, A2B, FIR, fft_size);
     end
   end
   if exist("FIR")
-    if (size(FIR, 1) <= 6)
+    if (size(FIR, 1) <= 8)
       %% filter the B format signal through the common FIR filters
       for m = (1:measurements)
 	for b = 1:size(FIR, 1)
@@ -225,7 +277,7 @@ function [SP, FC, FE] = sphear_signal_power_all(SIG, FS, fmin, fmax, frefmin, fr
   ref = sphear_signal_power(SIG, FS, frefmin, frefmax, fft_size);
   norm = mean(ref(:,:));
 
-  %% we only normalize with respect to WXYZ, UV in the frequency range
+  %% we only normalize with respect to WXYZ, STUV in the frequency range
   %% where WXYZ are well behaved has low power (high pass filter due
   %% to the difference microphones) so we don't want to average that.
   %% If we do that the A2B coefficients are much higher than the ones
@@ -233,8 +285,8 @@ function [SP, FC, FE] = sphear_signal_power_all(SIG, FS, fmin, fmax, frefmin, fr
   %% then attenuation in the filtering stage after that)
   WXYZnorm = mean(ref(:,[1,2,3,4]));
   XYavg = mean(WXYZnorm([2,3]));
-  %% use XY average for UV
-  norm = cat(1, WXYZnorm', [XYavg], [XYavg])';
+  %% use XY average for STUV
+  norm = cat(1, WXYZnorm', [XYavg], [XYavg], [XYavg], [XYavg])';
 
   for i = (1:size(FC))
     %% (measurement,signal,band,power)
@@ -252,7 +304,7 @@ end
 %% 2 -> use cardinal measurements
 %% 3 -> use diagonal measurements
 
-function [W X Y Z U V] = sphear_bf_extract_filter_shapes(B, measurements)
+function [W X Y Z S T U V] = sphear_bf_extract_filter_shapes(B, measurements)
   %% cardinal and diagonal directions for four capsule array
   TETRA_CARDINALS = [1, 5, 9, 13];
   TETRA_DIAGONALS = [3, 7, 11, 15];
@@ -262,9 +314,14 @@ function [W X Y Z U V] = sphear_bf_extract_filter_shapes(B, measurements)
   %% best directions for octa WXY
   OCTA_FLAT = OCTA_DIAGONALS;
   OCTA_BRIGHT = OCTA_CARDINALS;
+  %% best directions for octa Z
+  OCTA_FLAT_Z = [OCTA_FLAT + 16, OCTA_FLAT + 32];
   %% best directions for octa UV
   U_DIRS = [1, 5, 9, 13];
   V_DIRS = [3, 7, 11, 15];
+  %% best directions for octa ST (up and down layers)
+  S_DIRS = [U_DIRS + 16, U_DIRS + 32];
+  T_DIRS = [V_DIRS + 16, V_DIRS + 32];
   %% all directions
   ALL_DIRS = [1:16];
 
@@ -273,10 +330,19 @@ function [W X Y Z U V] = sphear_bf_extract_filter_shapes(B, measurements)
     if (size(B, 2) > 4)
       %% 8 capsules
       %% averaging all measurements for W seems to do better
+      %%W_m = [OCTA_FLAT, OCTA_FLAT + 16, OCTA_FLAT + 32];
+      %%Y_m = [OCTA_FLAT, OCTA_FLAT + 16, OCTA_FLAT + 32];
+      %%X_m = [OCTA_FLAT, OCTA_FLAT + 16, OCTA_FLAT + 32];
       W_m = OCTA_FLAT;
       Y_m = OCTA_FLAT;
       X_m = OCTA_FLAT;
-      Z_m = W_m;
+      if (size(B, 1) > 16)
+	Z_m = OCTA_FLAT_Z;
+      else
+	Z_m = W_m;
+      end
+      S_m = S_DIRS;
+      T_m = T_DIRS;
       U_m = U_DIRS;
       V_m = V_DIRS;
     else
@@ -323,8 +389,17 @@ function [W X Y Z U V] = sphear_bf_extract_filter_shapes(B, measurements)
   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);
+    if (size(B, 1) > 16)
+      %% full 3d measurements
+      S = 1.0 ./ reshape(mean(B(S_m,5,:)), size(B(:,:,:))(3), 1);
+      T = 1.0 ./ reshape(mean(B(T_m,6,:)), size(B(:,:,:))(3), 1);
+    else
+      %% flat 2d measurements, no S/T filters (ie: flat 1.0 filter)
+      S = 1.0 ./ ones(size(B, 3), 1);
+      T = 1.0 ./ ones(size(B, 3), 1);
+    end
+    U = 1.0 ./ reshape(mean(B(U_m,7,:)), size(B(:,:,:))(3), 1);
+    V = 1.0 ./ reshape(mean(B(V_m,8,:)), size(B(:,:,:))(3), 1);
   end
 end
 
@@ -332,31 +407,72 @@ 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)
+function [A2B XY_gain, S_gain, T_gain, U_gain, V_gain] = 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
+    do_Z = true;
+  end
   %% (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, peak of first order lobes
   XY_avg = mean([mean(B_power([1,9], 2)), mean(B_power([5, 13], 3))]);
+  XY_gain = XY_avg/W_avg;
+  if (do_Z)
+    %% Z, +-elev weighted by elevation angle
+    %% FIXME: we just hardwire the elevation angle for now
+    Z_avg = mean(B_power([17:32,33:48],4)) / sin(43/180*pi);
+    Z_gain = Z_avg/W_avg;
+  else
+    Z_gain = 0.0;
+  end
   if (size(A2B, 2) > 4)
     %% UV, peak of second order lobes
-    UV_avg = mean([mean(B_power([1,5,9,13], 5)), mean(B_power([3,7,11,15], 6))]);
+    U_avg = mean(mean(B_power([1,5,9,13], 7)));
+    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)
+      %% 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;
+      T_avg = mean([mean(B_power([5,13] + 16, 6)), mean(B_power([5, 13] + 32, 6))]);
+      T_gain = (T_avg/W_avg) * s_t_val;
+    else
+      S_gain = 0.0;
+      T_gain = 0.0;
+    end
+  else
+    U_gain = 0.0;
+    V_gain = 0.0;
+  end
+  %% normalize the A2B matrix
+  if (do_Z)
+    A2B(:,[2,3]) = A2B(:,[2,3]) / XY_gain;
+    A2B(:,[4]) = A2B(:,[4]) / Z_gain;
+  else
+    A2B(:,[2,3,4]) = A2B(:,[2,3,4]) / XY_gain;
   end
-  %% 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, 2) > 4)
-    UV_gain = UV_avg/W_avg;
-    A2B(:,[5,6]) = A2B(:,[5,6]) / UV_gain;
+    if (mean(mean(A2B(:,[5,6]))) > 0)
+      A2B(:,[5]) = A2B(:,[5]) / S_gain;
+      A2B(:,[6]) = A2B(:,[6]) / T_gain;
+    end
+    A2B(:,[7]) = A2B(:,[7]) / U_gain;
+    A2B(:,[8]) = A2B(:,[8]) / V_gain;
   else
-    UV_gain = 1.0;
   end
   %% we want the W contribution of each capsule to be 1/8,
   %% so scale the whole matrix relative to mean W
   A2B(:,:) = A2B(:,:) / (mean(A2B(:, 1)) * 8);
   %% print gain scaling coefficients
-  printf("XY gain: %f, UV gain: %f\n", XY_gain, UV_gain);
+  printf("XY gain: %f, Z gain: %f, ST: %f|%f, UV: %f|%f\n", XY_gain, Z_gain, S_gain, T_gain, U_gain, V_gain);
 end
 
 %%
@@ -477,3 +593,174 @@ function sphear_plot_all_filters(FC, UC, fmin, fmax)
   hold off;
 end
 
+%%%
+%%% Plot B format signal for three elevation measurement angles
+%%%
+
+function sphear_plot_measurements_3d(SIG, FC, FE, FS, N, freq)
+  figure;
+  hold on;
+  axis([1, 16]);
+  grid on;
+  ca=get (gcf, "currentaxes");
+  set(ca,"fontweight","bold","linewidth",1.5)
+  %% find the approximate desired frequency band
+  for f = [1:size(FC)]
+    if (FC(f) <= freq)
+      band = f;
+    end
+  end
+  title(["B format power in frequency band (" sprintf("%.0f", FE(band)) " to " sprintf("%.0f", FE(band+1)) " Hz)"]);
+  %% and plot the powers
+  for i = [1,17,33]
+    if (i < 17)
+      wstyle = "ko-;W0;";
+      color = "r";
+      Xname = "X 0elev";
+      Yname = "Y 0elev";
+    else
+      if (i < 33)
+	wstyle = "ko--;WU;";
+	color = "b";
+	Xname = "X +elev";
+	Yname = "Y +elev";
+      else
+	wstyle = "ko:;WD;";
+	color = "g";
+	Xname = "X -elev";
+	Yname = "Y -elev";
+      end
+    end
+    %% W
+    plot(SIG((i:i+15),1,band), wstyle, "linewidth", 1.5);
+    %% X/Y
+    plot(SIG((i:i+15),[2],band), "color", color, sprintf("o-;%s;", Xname), "linewidth", 1.5);
+    plot(SIG((i:i+15),[3],band), "color", color, sprintf("o-;%s;", Yname), "linewidth", 1.5);
+    %% Z
+    pZ = plot(SIG((i:i+15),4,band), "co-;Z;", "linewidth", 1.5);
+  end
+  legend('location', 'southwestoutside');
+  hold off;
+end
+
+%%%
+%%% create X matrix for given elevation angles
+%%%
+
+function [M x y z s t u v az el] = sphear_X_matrix(n, up, down)
+  eps = 0.01;
+  %% create arrays of angles
+  az = (0:360/n:360-eps)';
+  %% map zero elevation into spherical harmonic components
+  el = 0;
+  for i = 1:size(az)
+    x(i) = cos(az(i)/180*pi) * cos(el/180*pi);
+    y(i) = sin(az(i)/180*pi) * cos(el/180*pi);
+    z(i) = sin(el/180*pi);
+    s(i) = cos(az(i)/180*pi) * sin((2*el)/180*pi);
+    t(i) = sin(az(i)/180*pi) * sin((2*el)/180*pi);
+    u(i) = cos((2*az(i))/180*pi) * cos(el/180*pi) * cos(el/180*pi);
+    v(i) = sin((2*az(i))/180*pi) * cos(el/180*pi) * cos(el/180*pi);
+  end
+  %% map positive elevation measurements
+  el = up;
+  for i = 1:size(az)
+    x(i + 16) = cos(az(i)/180*pi) * cos(el/180*pi);
+    y(i + 16) = sin(az(i)/180*pi) * cos(el/180*pi);
+    z(i + 16) = sin(el/180*pi);
+    s(i + 16) = cos(az(i)/180*pi) * sin((2*el)/180*pi);
+    t(i + 16) = sin(az(i)/180*pi) * sin((2*el)/180*pi);
+    u(i + 16) = cos((2*az(i))/180*pi) * cos(el/180*pi) * cos(el/180*pi);
+    v(i + 16) = sin((2*az(i))/180*pi) * cos(el/180*pi) * cos(el/180*pi);
+  end
+  %% map negative elevation measurements
+  el = down;
+  for i = 1:size(az)
+    x(i + 32) = cos(az(i)/180*pi) * cos(el/180*pi);
+    y(i + 32) = sin(az(i)/180*pi) * cos(el/180*pi);
+    z(i + 32) = sin(el/180*pi);
+    s(i + 32) = cos(az(i)/180*pi) * sin((2*el)/180*pi);
+    t(i + 32) = sin(az(i)/180*pi) * sin((2*el)/180*pi);
+    u(i + 32) = cos((2*az(i))/180*pi) * cos(el/180*pi) * cos(el/180*pi);
+    v(i + 32) = sin((2*az(i))/180*pi) * cos(el/180*pi) * cos(el/180*pi);
+  end
+  %% return big matrix of projections
+  M = [x' y' z' s' t' u' v'];
+end
+
+%%%
+%%% Fractional delay
+%%%
+%%% based on:
+%%% http://www.labbookpages.co.uk/audio/beamforming/fractionalDelay.html
+%%%
+%%% Also:
+%%% http://users.spa.aalto.fi/vpv/publications/vesan_vaitos/ch3_pt1_fir.pdf
+%%%
+%%% N = length of sinc function (should be odd)
+
+function [SINC] = sphear_delayed_sinc(N, delay)
+  center = floor(N / 2);
+  for t = (1:N)
+    x = (t - 1) - delay;
+    window = 0.54 - 0.46 * cos(2.0 * pi * (x + 0.5) / N);
+    SINC(t) = sinc(x - center) * window;
+  end
+end
+
+function [DSIG FIR] = sphear_delay_signal(SIG, delay, N, fft_size)
+  FIR = sphear_delayed_sinc(N, delay);
+  for m = [1:size(SIG, 1)]
+    for c = [1:size(SIG, 2)]
+      %% we add trailing zeros to let the filter ring
+      DSIG(m, c, :) = fftfilt(FIR, cat(1, reshape(SIG(m, c, :), 1, [])', zeros(2*size(FIR, 2), 1))', fft_size)(round(N/2):end);
+    end
+  end
+end
+
+%%%
+%%% Fix wobble delays in a set of measurements
+%%%
+
+function [FIXED] = sphear_fix_wobble(SIG, DELAY);
+  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);
+  end
+end
+
+%%%
+%%% Do a basic encoder design, static A2B matrix plus FIR filters
+%%%
+
+function [BFF A2B FIR W X Y Z U V] = sphear_calculate_design(fref_min, fref_max)
+  fft_size = 8192;
+  N_oct = 16;
+  %% read in data
+  [AF, FS] = sphear_af_read((1:16));
+  %% power in reference band
+  [R] = sphear_signal_power(AF, FS, fref_min, fref_max, fft_size);
+  %% static A2B matrix
+  [A2B, PV, COND] = sphear_calculate_a2b(R);
+  %% generate a B-format signal
+  [BF] = sphear_af_to_bf(AF, A2B);
+  %% calculate power of B-format signals
+  [B, FC, FE] = sphear_signal_power_all(BF, FS, 300, 20000, fref_min, fref_max, N_oct, fft_size);
+  %% extract the shapes of wxyz[uv] filters
+  %% average response in cardinal directions
+  [W X Y Z U V] = sphear_bf_extract_filter_shapes(B, 0);
+  %% design FIR filters
+  [FIR IR MPS] = sphear_design_fir_filters(FC, 512, FS, W, X, Y, Z, U, V, 400, 18, -10);
+  %% re-calculate the full B-format signal with the filters
+  [BF, BFF] = sphear_af_to_bf(AF, A2B, FIR, fft_size);
+  %% get the average power in the directions of the lobes for W, XY and UV
+  %% and normalize the A2B matrix so that ZY and UV have the same peak amplitude
+  %% as W
+  [A2B XY_gain UV_gain] = sphear_bf_normalize_gain(A2B, BFF, FS, fref_min, fref_max, fft_size);
+  %% re-calculate the B-format signal
+  [BF, BFF] = sphear_af_to_bf(AF, A2B, FIR, fft_size);
+end
+

calibration/src/sphear_design_fir_filters.m

@@ -44,7 +44,7 @@
 %% NOTES: using fir2 for now, firls calculates non-smooth filters with worse
 %% frequency response, remez does not converge in my tests
 
-function [FIR IR MPS] = sphear_design_fir_filters(F, n, SR, W, X, Y, Z, U, V, UV_lf_freq, UV_lf_slope, DC_gain_offset_db)
+function [FIR IR MPS] = sphear_design_fir_filters(F, n, SR, W, X, Y, Z, S, T, U, V, UV_lf_freq, UV_lf_slope, DC_gain_offset_db)
   if ! exist("Z")
     %% if Z was not supplied, just average X and Y data
     Z = (X + Y) / 2;
@@ -82,8 +82,10 @@ function [FIR IR MPS] = sphear_design_fir_filters(F, n, SR, W, X, Y, Z, U, V, UV
 	IR(4+size(U,2)+c,:) = fir2(n, freqs, [V(1,c), V(:,c)', V(end,c)]' .* UV_lpf);
       end
     else
-      IR(5,:) = fir2(n, freqs, [U(1), U', U(end)]' .* UV_lpf);
-      IR(6,:) = fir2(n, freqs, [V(1), V', V(end)]' .* UV_lpf);
+      IR(5,:) = fir2(n, freqs, [S(1), S', S(end)]' .* UV_lpf);
+      IR(6,:) = fir2(n, freqs, [T(1), T', T(end)]' .* UV_lpf);
+      IR(7,:) = fir2(n, freqs, [U(1), U', U(end)]' .* UV_lpf);
+      IR(8,:) = fir2(n, freqs, [V(1), V', V(end)]' .* UV_lpf);
     end
   end
 

calibration/src/sphear_octa_af_calibrate.m

@@ -29,7 +29,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 [] = sphear_octa_af_calibrate(measurements, keep, fade);
+function [] = sphear_octa_af_calibrate(measurements, keep, fade, suffix);
+  if (! exist("suffix"))
+    suffix = "";
+  end
   %% the measurement from the reference microphone is processed through
   %% DRC (Digital Room Correction) to create a convolution kernel that is
   %% used to equalize the measurements
@@ -64,10 +67,10 @@ function [] = sphear_octa_af_calibrate(measurements, keep, fade);
     printf("processing measurement %d\n", m);
     %% read in capsule signals and normalize them
     for c = 1:8
-      sig = wavread(["../data/current/af/a-format-" sprintf("%03d", m) "-" char(octanames(c)) ".wav"]);
+      sig = wavread(["../data/current/af/a-format-" sprintf("%s%03d", suffix, m) "-" char(octanames(c)) ".wav"]);
       cal_sig(:,c) = conv(trunc, sig);
-      wavwrite(cal_sig(:,c), Fs, ["../data/current/af/a-format-cal-" sprintf("%03d", m) "-" char(octanames(c)) ".wav"]);
+      wavwrite(cal_sig(:,c), Fs, ["../data/current/af/a-format-cal-" sprintf("%s%03d", suffix, m) "-" char(octanames(c)) ".wav"]);
     end
-    wavwrite(cal_sig, Fs, ["../data/current/af/a-format-cal-" sprintf("%03d", m) ".wav"]);
+    wavwrite(cal_sig, Fs, ["../data/current/af/a-format-cal-" sprintf("%s%03d", suffix, m) ".wav"]);
   end
 end