calibration common code: changes names of files to account for new Aliki, add some support for experimenting with array of filters for U/V instead of single filters like for WXYZ

calibration/src/sphear_calibration_common.m

@@ -40,9 +40,9 @@ function [AF, FS] = sphear_af_read(measurements, raw = false);
   else
     type = "cal-";
   end
-  if (exist([DATA CURRENT "af/a-format-cal-001-C5.wav"]))
+  if (exist([DATA CURRENT "af/a-format-cal-001-5.wav"]))
     %% assume it is an eight capsule microphone
-    capnames = cellstr(['C1'; 'C2'; 'C3'; 'C4'; 'C5'; 'C6'; 'C7'; 'C8']);
+    capnames = cellstr(['1'; '2'; '3'; '4'; '5'; '6'; '7'; '8']);
     capsules = 8;
   else
     %% otherwise a tetrahedral microphone
@@ -135,12 +135,68 @@ function [BF, BFF] = sphear_af_to_bf(AF, A2B, FIR, fft_size);
     end
   end
   if exist("FIR")
-    %% filter the B format signal through the four FIR filters
-    for m = (1:measurements)
-      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);
+    if (size(FIR, 1) <= 6)
+      %% filter the B format signal through the common FIR filters
+      for m = (1:measurements)
+	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
+      %% first four components use common filters, U/V use filter arrays
+      for m = (1:measurements)
+	for b = 1:4
+	  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
+      capsules = 8;
+      for m = (1:measurements)
+	for b = 5:6
+	  for c = 1:capsules
+	    fi = 4 + (b - 5)*capsules + c;
+	    if (c == 1)
+	      BFF(m,b,:) = fftfilt(reshape(FIR(fi,:), 1, [])', cat(1, reshape(BF(m,b,:), 1, [])', zeros(2*size(FIR(1,:))(2),1))', fft_size);
+	    else
+	      BFF(m,b,:) = reshape(BFF(m,b,:), size(BFF, 3), 1)' .+ fftfilt(reshape(FIR(fi,:), 1, [])', cat(1, reshape(BF(m,b,:), 1, [])', zeros(2*size(FIR(1,:))(2),1))', fft_size);
+	    end
+	  end
+	end
+      end
+    end
+  else
+    BFF = BF;
+  end
+end
+
+%% Render a calibrated capsule impulse response through UV matrix filters
+%%
+%% m: number of measurement to render
+%% A2B: A to B matrix
+%% FIR: optional FIR filters
+%% fft_size: FFT size for filtering
+
+function [UV, UVF] = sphear_af_to_bf_new(AF, A2B, FIR, fft_size);
+  measurements = size(AF)(1);
+  BF = zeros(measurements, size(A2B, 2), size(AF, 3));
+  for m = (1:measurements)
+    for c = (1:size(AF,2))
+      bfc(1,:,:) = reshape(AF(m,c,:), 1, []) .* reshape(A2B(c,:), 1, [])';
+      BF(m,:,:) = BF(m,:,:) .+ bfc(1,:,:);
+    end
+  end
+  if exist("FIR")
+    if (size(FIR, 1) <= 6)
+      %% filter the B format signal through four or six FIR filters
+      for m = (1:measurements)
+	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
+      %% filter the B format signal through single filters for first order
+      %% and array of filters for second order
     end
   else
     BFF = BF;
@@ -151,7 +207,7 @@ end
 %% Calculate B format power in logarithmic frequency bands
 %%
 
-function [SP, FC, FE] = sphear_signal_power_all(SIG, FS, fmin, fmax, frefmin, frefmax, Noct, fft_size, Fs, plt)
+function [SP, FC, FE] = sphear_signal_power_all(SIG, FS, fmin, fmax, frefmin, frefmax, Noct, fft_size)
   %% calculate number of frequency bands spaced Noct
   n = floor(log(fmax/fmin)/log(2^(1/Noct)));
   %% adjust width of band so they fit exactly in frequency range
@@ -299,4 +355,125 @@ function [A2B XY_gain, UV_gain] = sphear_bf_normalize_gain(A2B, B, FS, fmin, fma
   %% 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);
+end
+
+%%
+%% Calculate A2B coefficients in logarithmic frequency bands
+%%
+
+function [A2BF] = sphear_calculate_a2b_all(SIG, FS, fmin, fmax, frefmin, frefmax, Noct, fft_size)
+  %% calculate number of frequency bands spaced Noct
+  n = floor(log(fmax/fmin)/log(2^(1/Noct)));
+  %% adjust width of band so they fit exactly in frequency range
+  N = log(2)/(log(fmax/fmin)/n);
+  printf("%d frequency bands, %f bands per octave\n", n, N);
+
+  %% first center frequency
+  fcmin = fmin * 2^(1/(N * 2));
+  %% all center and edge frequencies
+  FC = fcmin * 2.^((0:n - 1)/N)';
+  FE = fmin * 2.^((0:n)/N)';
+
+  band = 1;
+  for i = (1:size(FC))
+    %% get power in band
+    R = sphear_signal_power(SIG, FS, FE(i), FE(i + 1), fft_size);
+    %% calculate a2b matrix in band
+    [a2b, pv, cond] = sphear_calculate_a2b(R);
+    A2BF(band,:,:) = a2b;
+    band = band + 1;
+  end
+end
+
+%%
+%% Plot coefficients of A2B matrices as a function of frequency
+%%
+
+function sphear_plot_a2b_all(FC, A2B, fmin, fmax)
+  hold on;
+  grid on;
+  axis([fmin, fmax]);
+  xticks = [300, 500, 1000, 2000, 3000, 5000, 7000, 10000, 15000, 20000];
+  set(gca, 'XTick', xticks);
+  set(gca, 'XTickLabel', sprintf('%d|', xticks));
+  ca=get (gcf, "currentaxes");
+  set(ca,"fontweight","bold","linewidth",1.5)
+  cmap = jet(size(A2B, 2));
+  %% find the max frequency bin
+  ifmin = 0;
+  ifmax = 0;
+  for i = (1:size(A2B, 1))
+    if (FC(i) <= fmax)
+      ifmax = i;
+    end
+  end
+  ifmax
+  FC(ifmax)
+  size(A2B, 1)
+  %% plot
+  semilogx(FC([1:ifmax]), 10*log10(abs(A2B([1:ifmax],:,:))), ".-", "linewidth", 1.5);
+  l = legend("C1", "C2", "C3", "C4", "C5", "C6", "C7", "C8");
+  %%legend(l, 'location', 'west', 'outside');
+  hold off;
+end
+
+%%%
+%%% Merge 2nd order low frequency filter shape with overall filter shape
+%%%
+
+function [UC, VC] = sphear_merge_filters(FC, FE, U, V, A2BF, fmax)
+  %% find max frequency index
+  imax = 0;
+  for f = 1:size(FC, 1)
+    if (FC(f) < fmax)
+      imax = f;
+    end
+  end
+  %% merge filters below the maximum frequency
+  for f = 1:size(FC, 1)
+    if (FC(f) < fmax)
+      for c = 1:size(A2BF, 2)
+	%%UC(f,c,:) = U(f) + real(10 * log10(abs(A2BF(f,c,5) / A2BF(imax,c,5))));
+	%%VC(f,c,:) = V(f) + real(10 * log10(abs(A2BF(f,c,6) / A2BF(imax,c,6))));
+	UC(f,c,:) = U(f) + real(10 * log10(abs(A2BF(imax,c,5) / A2BF(f,c,5))));
+	VC(f,c,:) = V(f) + real(10 * log10(A2BF(imax,c,6) / abs(A2BF(f,c,6))));
+      end
+    else
+      for c = 1:size(A2BF, 2)
+	UC(f,c,:) = U(f);
+	VC(f,c,:) = V(f);
+      end
+    end
+  end
 end
+
+%%%
+%%% Plot component filters
+%%%
+
+function sphear_plot_all_filters(FC, UC, fmin, fmax)
+  hold on;
+  grid on;
+  axis([fmin, fmax]);
+  xticks = [300, 500, 1000, 2000, 3000, 5000, 7000, 10000, 15000, 20000];
+  set(gca, 'XTick', xticks);
+  set(gca, 'XTickLabel', sprintf('%d|', xticks));
+  ca=get (gcf, "currentaxes");
+  set(ca,"fontweight","bold","linewidth",1.5)
+  %% find the max frequency bin
+  ifmin = 0;
+  ifmax = 0;
+  for i = (1:size(UC, 1))
+    if (FC(i) <= fmax)
+      ifmax = i;
+    end
+  end
+  %% plot
+  semilogx(FC([1:ifmax]), UC([1:ifmax],:,:), ".-", "linewidth", 1.5);
+  l = legend("C1", "C2", "C3", "C4", "C5", "C6", "C7", "C8");
+  %%legend(l, 'location', 'west', 'outside');
+  hold off;
+end
+