add common short support functions for calibration, load before starting the process

calibration/src/sphear_calibration_common.m

@@ -0,0 +1,174 @@
+%% -*- mode: octave -*-
+%%
+%% *SpHEAR Project
+%%
+%% Common functions for A2B encoder design
+%%
+%% Copyright 2016-2017, Fernando Lopez-Lezcano, All Rights Reserved
+%% nando@ccrma.stanford.edu
+%%
+%%  This program is free software: you can redistribute it and/or modify
+%%  it under the terms of the GNU General Public License as published by
+%%  the Free Software Foundation, either version 3 of the License, or
+%%  (at your option) any later version.
+%%
+%%  This program is distributed in the hope that it will be useful,
+%%  but WITHOUT ANY WARRANTY; without even the implied warranty of
+%%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+%%  GNU General Public License for more details.
+%%
+%%  You should have received a copy of the GNU General Public License
+%%  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+%% make sure octave treats this as a script file
+1;
+
+pkg load signal;
+
+%%
+%% Read an A-Format signal for specified measurements
+%%
+
+global DATA = "../data/";
+global CURRENT = "current/";
+
+function [AF, FS] = sphear_af_read(measurements, raw = false);
+  global DATA;
+  global CURRENT;
+  if (raw)
+    type = "";
+  else
+    type = "cal-";
+  end
+  if (exist([DATA CURRENT "af/a-format-cal-001-C5.wav"]))
+    %% assume it is an eight capsule microphone
+    capnames = cellstr(['C1'; 'C2'; 'C3'; 'C4'; 'C5'; 'C6'; 'C7'; 'C8']);
+    capsules = 8;
+  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
+      [A, FS] = wavread([DATA CURRENT "af/a-format-" type sprintf("%03d", m) "-" char(capnames(c)) ".wav"]);
+      AF(m, c, :) = A;
+    end
+  end
+end
+
+%%
+%% Calculate power for signal in a frequency range
+%%
+%% (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)
+  %% 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)
+      if (s == 4 && cardioid)
+	spectrum = fft((SIG(m, 1, :) + 1.0 * SIG(m, 2, :)) / 2, fft_size);
+      else
+	spectrum = fft(SIG(m, s, :), fft_size);
+      end
+      %% 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);
+  %% normalize matrix
+  A2B = A2B ./ mean(mean(abs(A2B)));
+  %% return condition number
+  COND = cond(A2B);
+end
+
+%% Render a calibrated capsule impulse response through an A2B matrix
+%%
+%% m: number of measurement to render
+%% A2B: A to B matrix
+%% FIR: optional FIR filters
+%% fft_size: FFT size for filtering
+
+function [BF, BFF] = sphear_af_to_bf(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")
+    %% filter the B format signal through the four FIR filters
+    for m = (1:measurements)
+      for b = 1:4
+	%% 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
+  end
+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)
+  %% 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)';
+
+  %% normalize power with respect to reference frequency band
+  %% (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;
+  end
+end
+
+%%
+%% Extract average shapes of B format filters
+%%
+
+%% assuming 16 measurements equally spaced around the microphone
+TETRA_CARDINALS = [1, 5, 9, 13];
+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)
+  %% extract average response shape for W
+  W = 1.0 ./ reshape(mean(B(measurements,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);
+end
+