calibration: assume all impulse responses have been trimmed and windowed properly, do not try to determine impulse location; many changes in parameters and parameter order; do a running average of filter coefficients when designing the filter transition; determine the transition frequency automatically based on the maximum condition number of the A2B matrices

calibration/src/sphear_a2b_all.m

@@ -4,19 +4,18 @@
 %% sphear_a2b_all: calculate A2B matrixes for frequency ranges
 %%
 %% usage:
-%%   [F FE A2BF Fs] = sphear_a2b_all(m, fmin, fmax, Noct, Nsamp, Npre, plot_on)
+%%   [F FE Fs A2BF PVF COND TF] = sphear_a2b_all(m, fmin, fmax, Noct, Nfft, plot_on)
 %%
 %% m: vector of measurement numbers to take into account
 %% fmin, fmax: lower and upper frequency limits
 %% Noct: width of frequency ranges in fractions of an octave, this
 %%       will be changed slightly so that the frequency range is
 %%       covered by equally spaced frequency bands
-%% Nsamp: total number of time domain samples to use
-%% Npre: samples to use before the main impulse peak
+%% Nfft: fft size for calculating power
 %%
 %% optional parameters:
 %%
-%% dplot: plot polar response for each frequency band (default is 0)
+%% plot_on: plot polar response for each frequency band (default is 0)
 %%
 %% Copyright 2016, Fernando Lopez-Lezcano, All Rights Reserved
 %% nando@ccrma.stanford.edu
@@ -34,7 +33,7 @@
 %%  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 [F FE Fs A2BF PVF] = sphear_a2b_all(m, fmin, fmax, Noct, Nsamp, Npre, plot_on)
+function [F FE Fs A2BF PVF COND TF] = sphear_a2b_all(m, fmin, fmax, Noct, Nfft, plot_on)
 
   %% do not plot by default
   if (!exist("plot_on"))
@@ -54,21 +53,28 @@ function [F FE Fs A2BF PVF] = sphear_a2b_all(m, fmin, fmax, Noct, Nsamp, Npre, p
   %% all center and edge frequencies
   fc = fcmin * 2.^((0:n - 1)/N)';
   fe = fmin * 2.^((0:n)/N)';
-  
+
   for i = (1:size(fc))
-    [R srate] = sphear_af_power(m, fe(i), fe(i + 1), Nsamp, Npre, 1.1);
+    [R srate] = sphear_af_power(m, fe(i), fe(i + 1), Nfft);
     [A2B pv max_err] = A2B_matrix(R);
     A2BF(i,:,:) = A2B;
     PVF(i,:,:) = pv;
     BF = R * A2B;
+    COND(i) = cond(A2B);
 
-    printf("%2d: %5d -> [%5d %5d] A2B cond=%f\n", i, fc(i), fe(i), fe(i + 1), cond(A2B));
+    printf("%2d: %5d -> [%5d %5d] A2B cond=%f\n", i, fc(i), fe(i), fe(i + 1), COND(i));
 
     if plot_on
       sphear_bf_plot(BF, fe(i), fe(i + 1), 1.1);
     end
   end
 
+  %% calculate the transition frequency of the array
+  [value, index] = max(COND);
+  TF = fc(index) / 2;
+  printf("max condition %fHz, cf %f, fe %f, transition %f\n", value, fc(index), fe(index), fc(index)/2);
+  
+  
   F = fc;
   FE = fe;
   Fs = srate;

calibration/src/sphear_a2b_fir_all.m

@@ -4,16 +4,16 @@
 %% sphear_a2b_all: calculate A2B matrixes for frequency ranges
 %%
 %% usage:
-%%   [F FE A2BF] = sphear_a2b_all(m, fmin, fmax, Noct, Fs, Nsamp, Npre, plt)
+%%   [F FE A2BF] = sphear_a2b_all(m, fmin, fmax, Noct, Nfft, Fs, plt)
 %%
 %% m: vector of measurement numbers to take into account
 %% fmin, fmax: lower and upper frequency limits
 %% Noct: width of frequency ranges in fractions of an octave, this
 %%       will be changed slightly so that the frequency range is
 %%       covered by equally spaced frequency bands
+%% Nfft: fft size
 %% Fs: sampling rate
-%% Nsamp: time domain samples to use
-%% Npre: samples to use before the main impulse peak
+%% plt: plot
 %%
 %% Copyright 2016, Fernando Lopez-Lezcano, All Rights Reserved
 %% nando@ccrma.stanford.edu

calibration/src/sphear_af_power.m

@@ -4,11 +4,10 @@
 %% sphear_af_power: return array of power for all capsules in a range of frequencies
 %%
 %% usage:
-%%   [R Fs] = sphear_af_power(fmin, fmax, N, Npre)
+%%   [R Fs] = sphear_af_power(fmin, fmax, Nfft)
 %%
 %% fmin, fmax: range of frequencies to consider
-%% Npre: number of samples before the main impulse
-%% N: total number of samples to analyze
+%% Nfft: fft size
 %%
 %% The names of reference microphone impulse response and capsule measurements are
 %% currently hardwired into this file, this needs to be fixed and turned into
@@ -30,20 +29,16 @@
 %%  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 [R Fs] = sphear_af_power(measurements, fmin, fmax, N, Npre)
+function [R Fs] = sphear_af_power(measurements, fmin, fmax, Nfft)
   pkg load signal;
 
-  %% read in reference microphone, select samples to process
+  %% read in reference microphone
   [ref, Fs] = wavread('../data/current/ref/reference-001-1.wav');
-  [peaks, locs] = findpeaks(abs(ref));
-  [~, index] = max(peaks);
-  start = locs(index) - Npre;
-  range = start:(start + N - 1);
-  ref_fft = fft(ref(range), N);
+  ref_fft = fft(ref, Nfft);
   
   %% indexes in array for frequency bins
-  imin = round((N/2) * fmin / (Fs/2));
-  imax = round((N/2) * fmax / (Fs/2));
+  imin = round((Nfft/2) * fmin / (Fs/2));
+  imax = round((Nfft/2) * fmax / (Fs/2));
 
   capnames = cellstr(['LF'; 'RF'; 'LB'; 'RB']);
 
@@ -58,7 +53,7 @@ function [R Fs] = sphear_af_power(measurements, fmin, fmax, N, Npre)
       %% anything outside is undefined
       %%
       capsule = wavread(["../data/current/af/a-format-0" sprintf("%02.0f", m) "-" char(capnames(c)) ".wav"]);
-      equalized = fft(capsule(range), N) ./ ref_fft;
+      equalized = fft(capsule, Nfft) ./ ref_fft;
       
       %% return total power of signal in frequency range
       R(i, c) = sqrt(sum(equalized(imin:imax).*conj(equalized(imin:imax)))/size(equalized(imin:imax),1));

calibration/src/sphear_bf_plot.m

@@ -54,7 +54,7 @@ function [BF] = sphear_bf_plot(BF, fmin, fmax, rmax)
     axis([-rmax, rmax, -rmax, rmax]);
     ca=get (gcf, "currentaxes");
     set(ca,"fontweight","bold","linewidth",1.5)
-    title([char(bname(b)) ": " sprintf("%05.0d", fmin) " to " sprintf("%05.0d", fmax) " Hz"]);
+    title([char(bname(b)) ": " sprintf("%d", round(fmin)) " to " sprintf("%d", round(fmax)) " Hz"]);
 
     if (strcmp(plot_on, "epsc") || strcmp(plot_on, "pngc"))
       %% color plot

calibration/src/sphear_common.m

@@ -31,14 +31,13 @@ pkg load signal;
 %% range: range of samples to process in file
 %% Nfft: size of FFT to use
 
-function [bf] = sphear_a_to_b(m, FIR, range, Nfft);
+function [bf] = sphear_a_to_b(m, FIR, Nfft);
   capnames = cellstr(['LF'; 'RF'; 'LB'; 'RB']);
   %% read in capsule signals and apply FIR filter
   for c = 1:4
     cs(c,:) = wavread(["../data/current/af/a-format-0" sprintf("%02.0f", m) "-" char(capnames(c)) ".wav"]);
-    cr(c,:) = tukeywin(size(range)(2), 0.0) .* cs(c,:)'(range);
     for b = 1:4
-      bs(c,b,:) = fftfilt(reshape(FIR(b,c,:), 1, [])', cr(c,:), Nfft);
+      bs(c,b,:) = fftfilt(reshape(FIR(b,c,:), 1, [])', cs(c,:), Nfft);
     end
   end
   %% add filtered capsule output to b format components

calibration/src/sphear_fir_design.m

@@ -63,12 +63,13 @@ function [FIR] = sphear_fir_design(F, A2BF, n, maxboost, Fs, Fl, plot_on, W, X,
   %% create a frequency response profile in dB
   fr_profile = interp1([F(1), F(fl_idx), F(end)], [0, 0, 0], F);
 
+  dmax = size(A2BF(:,1,1))(1);
   %% all b format components
   for b = 1:4
     %% all capsules
     for c = 1:4
       %% all frequency bands
-      for d = 1:size(A2BF(:,1,1))(1)
+      for d = 1:dmax
 	%% create an array of points to design the filter
 	if (F(d) <= Fl)
 	  %% below the limit frequency we use the A2B coefficients
@@ -77,10 +78,19 @@ function [FIR] = sphear_fir_design(F, A2BF, n, maxboost, Fs, Fl, plot_on, W, X,
 	  A2BF(d, c, b) = vlast * 10^(fr_profile(d)/20);
 	  dlast = d;
 	  if (exist("W") == 1)
-	    wlast = W(d);
-	    xlast = X(d);
-	    ylast = Y(d);
-	    zlast = Z(d);
+	    if (d == 1)
+	      wlast = W(d);
+	      xlast = X(d);
+	      ylast = Y(d);
+	      zlast = Z(d);
+	    else
+	      if (d > 1)
+		wlast = (wlast + W(d)) / 2;
+		xlast = (xlast + X(d)) / 2;
+		ylast = (ylast + Y(d)) / 2;
+		zlast = (zlast + Z(d)) / 2;
+	      end
+	    end
 	  end
 	else
 	  %% above the limit frequency we insert the W, X, Y and Z

calibration/src/sphear_fir_render.m

@@ -4,13 +4,12 @@
 %% sphear_fir_render: render impulse responses throught the FIR filter A2B matrix
 %%
 %% usage:
-%%   [] = sphear_fir_render(FIR, fmin, fmax, ymin, ymax, N, Npre, Nfft, plot_on, smooth)
+%%   [] = sphear_fir_render(FIR, fmin, fmax, ymin, ymax, Nfft, plot_on, smooth)
 %%
 %% FIR: A2B matrix of FIR filters
 %% fmin, fmax: range of frequencies to use
 %% ymin, ymax: y range to plot
-%% N: number of samples to filter
-%% Npre: number of samples before the main impulse peak
+%% Fs: sampling rate
 %% Nfft: size of fft used to calculate frequency response
 %%
 %% optional parameters:
@@ -35,7 +34,7 @@
 %%  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 [BFF] = sphear_fir_render(FIR, fmin, fmax, ymin, ymax, N, Npre, Nfft, plot_on, smooth)
+function [BFF] = sphear_fir_render(FIR, fmin, fmax, ymin, ymax, Fs, Nfft, plot_on, smooth)
 
   %% no smoothing or plotting by default
   if (!exist("smooth"))
@@ -49,19 +48,9 @@ function [BFF] = sphear_fir_render(FIR, fmin, fmax, ymin, ymax, N, Npre, Nfft, p
     end
   end
 
-  [ref, Fs] = wavread('../data/current/ref/reference-001-1.wav');
-  %% find highest peak location and define sample range
-  [index, value] = findmaxpeak(ref);
-  r_start = index - Npre;
-  r_end = r_start + N -1;
-  range = r_start:r_end;
-  ref_fft = fft(ref(range), N);
-  printf("reference: peak %f at %d, sample range [%d %d], sampling rate: %d\n",
-	 value, index, r_start, r_end, Fs);
-  
   for m = 1:16
     %% render b format components
-    bf(m,:,:) = sphear_a_to_b(m, FIR, range, N);
+    bf(m,:,:) = sphear_a_to_b(m, FIR, Nfft);
 
     %% calculate B format frequency responses
     for b = 1:4

calibration/src/sphear_fir_render_directions.m

@@ -3,12 +3,12 @@
 %% sphear_fir_render_directions: render impulse responses throught the FIR filter A2B matrix
 %%
 %% usage:
-%%   [] = sphear_fir_render_directions(FIR, fmin, fmax, N, Npre, Nfft, plot_on, smooth)
+%%   [] = sphear_fir_render_directions(FIR, fmin, fmax, ymin, ymax, Nfft, plot_on, smooth)
 %%
 %% FIR: A2B matrix of FIR filters
 %% fmin, fmax: range of frequencies to use
-%% N: number of samples to filter
-%% Npre: number of samples before the main impulse peak
+%% ymin, ymax: vertical scale for plot
+%% Fs: sampling rate
 %% Nfft: size of FFT
 %%
 %% optional parameters:
@@ -33,7 +33,7 @@
 %%  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 [BFF] = sphear_fir_render_directions(FIR, fmin, fmax, N, Npre, Nfft, plot_on, smooth)
+function [BFF] = sphear_fir_render_directions(FIR, fmin, fmax, ymin, ymax, Fs, Nfft, plot_on, smooth)
 
   %% no smoothing or plotting by default
   if (!exist("smooth"))
@@ -47,15 +47,6 @@ function [BFF] = sphear_fir_render_directions(FIR, fmin, fmax, N, Npre, Nfft, pl
     end
   end
 
-  [ref, Fs] = wavread('../data/current/ref/reference-001-1.wav');
-  %% find highest peak location and define sample range
-  [index, value] = findmaxpeak(ref);
-  r_start = index - Npre;
-  r_end = r_start + N -1;
-  range = r_start:r_end;
-  ref_fft = fft(tukeywin(size(range)(2), 0.02) .* ref(range), N);
-  printf("reference: peak %f at %d, sample range [%d %d]\n", value, index, r_start, r_end);
-  
   %% plot if requested
   if (plot_on != 0)
     if (strcmp(plot_on, "epsc") || strcmp(plot_on, "pngc"))
@@ -76,8 +67,8 @@ function [BFF] = sphear_fir_render_directions(FIR, fmin, fmax, N, Npre, Nfft, pl
     f = figure;
     hold on;
     grid on;
-    axis([fmin, fmax, -40, 0]);
-    tics('y', [6, 0, -3, -6, -9, -12, -15, -20, -25, -30]);
+    axis([fmin, fmax, ymin, ymax]);
+    tics('y', [6, 3, 0, -3, -6, -9, -12, -15, -20, -25, -30]);
     tics('x', [300, 500, 1000, 2000, 3000, 5000, 7000, 10000, 15000, 20000]);
     ca=get (gcf, "currentaxes");
     set(ca,"fontweight","bold","linewidth",1.5)
@@ -85,7 +76,7 @@ function [BFF] = sphear_fir_render_directions(FIR, fmin, fmax, N, Npre, Nfft, pl
     idx = 1;
     for m = [1, 5, 9, 13, 3, 7, 11, 15]
       %% render b format components
-      bf(idx,:,:) = sphear_a_to_b(m, FIR, range, N);
+      bf(idx,:,:) = sphear_a_to_b(m, FIR, Nfft);
 
       %% calculate B format frequency responses
       for b = 1:4