add a bunch of support functions for new encoder design process

calibration/src/sphear_af_plot_frequency_response.m

@@ -0,0 +1,81 @@
+%% -*- mode: octave -*-
+%%
+%% *SpHEAR Project
+%%
+%% sphear_af_plot_frequency_response: plot frequency response of B-Format signals
+%%
+%% usage:
+%%   [] = sphear_af_plot_frequency_response(SIG, fmin, fmax, ymin, ymax, Fs, Nfft, smooth)
+%%
+%% SIG: signals for each measurement
+%% fmin, fmax: range of frequencies to use
+%% ymin, ymax: y range to plot
+%% Fs: sampling rate
+%% Nfft: size of fft used to calculate frequency response
+%%
+%% optional parameters:
+%%
+%% smooth: number of samples to average, default is 0
+%%
+%% 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 <htp://www.gnu.org/licenses/>.
+
+function [] = sphear_af_plot_frequency_response(SIG, fmin, fmax, ymin, ymax, Fs, Nfft, smooth)
+
+  %% no smoothing or plotting by default
+  if (!exist("smooth"))
+    smooth = 0;
+  end
+  for m = 1:16
+    %% calculate frequency responses
+    for c = 1:size(SIG,2)
+      if (smooth > 0)
+	sfft(m,c,:) = fastsmooth(20*log10(abs(fft(SIG(m,c,:), Nfft))), smooth, 2, 1);
+      else
+	sfft(m,c,:) = 20*log10(abs(fft(SIG(m,c,:), Nfft)));
+      end
+    end
+    %% color plot
+    copts = ['k', 'r', 'b', 'g', ':g', ':b', ':r', ':k'];
+    cext = "_color";
+
+    f = figure;
+    hold on;
+    fr = (1:Nfft/2) * Fs/Nfft;
+    grid on;
+    angle = 360 / 16 * (m - 1);
+    %% we add the title as a subfigure in the paper
+    %% title(sprintf("azimuth = %.1f degrees", angle));
+    axis([fmin, fmax, ymin, ymax]);
+    tics('y', [6, 3, 0, -3, -6, -9, -12, -15, -20, -25, -30, -40, -50]);
+    tics('x', [300, 500, 1000, 2000, 3000, 5000, 7000, 10000, 15000, 20000]);
+    ca=get (gcf, "currentaxes");
+    set(ca,"fontweight","bold","linewidth",1.5)
+    pbaspect([16, 9, 1])
+    for c = 1:size(SIG,2)
+      semilogx(fr, sfft(m,c,1:end/2), copts(c), "linewidth", 1.5);
+    end
+    if (size(SIG,2) > 4)
+      %%l = legend('C1', 'C2', 'C3', 'C4', 'C5', 'C6', 'C7', 'C8');
+    else
+      %%l = legend('C1', 'C2', 'C3', 'C4');
+    end
+    %%legend(l, 'location', 'west');
+    hold off;
+    print(f, '-deps', '-tight', '-color', ["../data/current/figures/sphear_bf_azimuth_" sprintf("%03.0f", angle) cext ".eps"]);
+  end
+end
+

calibration/src/sphear_bf_plot_frequency_response.m

@@ -0,0 +1,110 @@
+%% -*- mode: octave -*-
+%%
+%% *SpHEAR Project
+%%
+%% sphear_bf_plot_frequency_response: plot frequency response of B-Format signals
+%%
+%% usage:
+%%   [] = sphear_bf_plot_frequency_response(BF, fmin, fmax, ymin, ymax, Fs, Nfft, plot_on, smooth)
+%%
+%% BF: B-Format signals for each measurement
+%% fmin, fmax: range of frequencies to use
+%% ymin, ymax: y range to plot
+%% Fs: sampling rate
+%% Nfft: size of fft used to calculate frequency response
+%%
+%% optional parameters:
+%%
+%% plot_on: plot frequency response of components, default is "epsc"
+%%          can be one of 0, "eps", "epsc" (color), "png", "pngc" (color)
+%% smooth: number of samples to average, default is 0
+%%
+%% 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/>.
+
+function [] = sphear_bf_plot_frequency_response(BF, fmin, fmax, ymin, ymax, Fs, Nfft, plot_on, smooth)
+
+  %% no smoothing or plotting by default
+  if (!exist("smooth"))
+    smooth = 0;
+  end
+  if (!exist("plot_on"))
+    plot_on = "epsc";
+  else
+    if (plot_on == 1)
+      plot_on = "epsc"
+    end
+  end
+
+  for m = 1:16
+    %% calculate B format frequency responses
+    for b = 1:4
+      if (smooth > 0)
+	bfft(m,b,:) = fastsmooth(20*log10(abs(fft(BF(m,b,:), Nfft))), smooth, 2, 1);
+      else
+	bfft(m,b,:) = 20*log10(abs(fft(BF(m,b,:), Nfft)));
+      end
+    end
+
+    %% plot if requested
+    if (plot_on != 0)
+      if (strcmp(plot_on, "epsc") || strcmp(plot_on, "pngc"))
+	%% color plot
+	wopts = 'k;W;';
+	xopts = 'r;X;';
+	yopts = 'b;Y;';
+	zopts = 'g;Z;';
+	copt = '-color';
+	cext = "_color";
+      else
+	%% mono plot
+	wopts = '-;W;';
+	xopts = '--;X;';
+	yopts = '-.;Y;';
+	zopts = ':;Z;';
+	copt = '-mono';
+	cext = "_mono";
+      end
+      f = figure;
+      hold on;
+      fr = (1:Nfft/2) * Fs/Nfft;
+      grid on;
+      angle = 360 / 16 * (m - 1);
+      %% we add the title as a subfigure in the paper
+      %% title(sprintf("azimuth = %.1f degrees", angle));
+      axis([fmin, fmax, ymin, ymax]);
+      tics('y', [6, 3, 0, -3, -6, -9, -12, -15, -20, -25, -30, -40, -50]);
+      tics('x', [300, 500, 1000, 2000, 3000, 5000, 7000, 10000, 15000, 20000]);
+      ca=get (gcf, "currentaxes");
+      set(ca,"fontweight","bold","linewidth",1.5)
+      pbaspect([16, 9, 1])
+      semilogx(fr, bfft(m,1,1:end/2), wopts, "linewidth", 1.5);
+      semilogx(fr, bfft(m,2,1:end/2), xopts, "linewidth", 1.5);
+      semilogx(fr, bfft(m,3,1:end/2), yopts, "linewidth", 1.5);
+      semilogx(fr, bfft(m,4,1:end/2), zopts, "linewidth", 1.5);
+      l = legend('W', 'X', 'Y', 'Z');
+      legend(l, 'location', 'west');
+      hold off;
+      if (strcmp(plot_on, 'eps') || strcmp(plot_on, 'epsc') || plot_on == 1)
+	print(f, '-deps', '-tight', copt, ["../data/current/figures/sphear_bf_azimuth_" sprintf("%03.0f", angle) cext ".eps"]);
+      else
+	print(f, '-dpng', copt, ["../data/current/figures/sphear_bf_azimuth_" sprintf("%03.0f", angle) cext ".png"]);
+      end
+    end
+  end
+
+end
+

calibration/src/sphear_bf_plot_polar_response.m

@@ -0,0 +1,76 @@
+%% -*- mode: octave -*-
+%%
+%% *SpHEAR Project
+%%
+%% sphear_bf_plot_polar_response: plot polar response of B format components
+%%
+%% usage:
+%%   sphear_bf_plot_polar_response(B, frefmin, frefmax, fmin, fmax, Noct, Nfft, Fs, cardioid = 0)
+%%
+%% B: B format signal
+%% frefmin, frefmax: reference frequencies for power reference
+%% 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
+%%
+%% optional:
+%%
+%% cardioid: 1 -> replace Z component by a virtual cardioid pointing forward
+%%
+%% 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/>.
+
+function sphear_bf_plot_polar_response(B, frefmin, frefmax, fmin, fmax, Noct, Nfft, Fs, cardioid = 0)
+
+  %% color plot
+  copt = '-color';
+  cext = "_color";
+
+  plotz = 0;
+  
+  %% 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 amplitude with respect to W between 600 and 1200 Hz
+  ref = sphear_signal_power(B, Fs, frefmin, frefmax, Nfft);
+  norm = mean(ref(:,1));
+  
+  f = figure;
+
+  idx = 1;
+  for i = (1:size(fc))
+    printf("fc = %d\n", fc(i));
+    BP = sphear_signal_power(B, Fs, fe(i), fe(i + 1), Nfft, cardioid) / norm;
+    sphear_bf_plot(BP, fe(i), fe(i + 1), 1.25, cardioid);
+    idx = idx + 1;
+  end
+
+  %% print(f, '-deps','-tight',copt,["figures/sphear_cardioid_polar.eps"]);
+
+end
+

calibration/src/sphear_capsule_plot_frequency_response.m

@@ -0,0 +1,89 @@
+%% -*- mode: octave -*-
+%%
+%% *SpHEAR Project
+%%
+%% sphear_capsule_plot_frequency_response: plot frequency response of capsule signals
+%%
+%% usage:
+%%   [] = sphear_capsule_plot_frequency_response(SIG, fmin, fmax, ymin, ymax, Nfft, plot_on, smooth)
+%%
+%% SIG: signals for each measurement
+%% fmin, fmax: range of frequencies to use
+%% ymin, ymax: y range to plot
+%% Fs: sampling rate
+%% Nfft: size of fft used to calculate frequency response
+%%
+%% optional parameters:
+%%
+%% plot_on: plot frequency response of components, default is "epsc"
+%%          can be one of 0, "eps", "epsc" (color), "png", "pngc" (color)
+%% smooth: number of samples to average, default is 0
+%%
+%% 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 <htp://www.gnu.org/licenses/>.
+
+function [] = sphear_capsule_plot_frequency_response(SIG, measurements, fmin, fmax, ymin, ymax, Fs, Nfft, smooth)
+
+  %% no smoothing or plotting by default
+  if (!exist("smooth"))
+    smooth = 0;
+  end
+  for c = 1:size(SIG,2)
+    %% calculate frequency responses
+    for m = measurements
+      if (smooth > 0)
+	sfft(m,c,:) = fastsmooth(20*log10(abs(fft(SIG(m,c,:), Nfft))), smooth, 2, 1);
+      else
+	sfft(m,c,:) = 20*log10(abs(fft(SIG(m,c,:), Nfft)));
+      end
+    end
+    %% color plot
+    copts = ['k', 'r', 'b', 'g', ':g', ':b', ':r', ':k'];
+    cext = "_color";
+
+    f = figure;
+    hold on;
+    fr = (1:Nfft/2) * Fs/Nfft;
+    grid on;
+    angle = 360 / 16 * (m - 1);
+    %% we add the title as a subfigure in the paper
+    %% title(sprintf("azimuth = %.1f degrees", angle));
+    axis([fmin, fmax, ymin, ymax]);
+    tics('y', [6, 3, 0, -3, -6, -9, -12, -15, -20, -25, -30, -40, -50]);
+    tics('x', [300, 500, 1000, 2000, 3000, 5000, 7000, 10000, 15000, 20000]);
+    ca=get (gcf, "currentaxes");
+    set(ca,"fontweight","bold","linewidth",1.5)
+    pbaspect([16, 9, 1])
+    for m = measurements
+      if (m == 1)
+	col = "k";
+      else
+	if (m < 8)
+	  col = "r";
+	else
+	  if (m == 9)
+	    col = "g";
+	  else
+	    col = "b";
+	  end
+	end
+      end
+      semilogx(fr, sfft(m,c,1:end/2), col, "linewidth", 1.5);
+    end
+    hold off;
+  end
+end
+

calibration/src/sphear_design_fir_filters.m

@@ -0,0 +1,95 @@
+%% -*- mode: octave -*-
+%%
+%% *SpHEAR Project
+%%
+%% sphear_design_fir_filters: design FIR filters based on W, X, Y and Z shapes
+%%
+%% usage:
+%%   [FIR] = sphear_design_fir_filters(F, n, maxboost, Fs, plot_on,
+%%                                     shelf_gain, shelf_freq, W, X, Y, Z)
+%%
+%% F: array of filter section center frequencies
+%% n: order of the FIR filter being created
+%% maxboost: max high and low frequency boost in dB of the final filter
+%% Fs: sampling rate
+%% plot_on: plot frequency response of components, default is "epsc"
+%%          can be one of 0, "eps", "epsc" (color), "png", "pngc" (color)
+%%
+%% shelf_gain: optional high pass shelf filter, gain in dB (0.0 dB turns it off)
+%% shelf_freq: corner frequency in Hertz
+%%
+%% for the first run do not supply W/X/Y/Z, render the filters and then run
+%% sphear_fir_extract to extract the W, X, Y and Z filter shapes, finally run
+%% this again with the extra parameters below
+%%
+%% W: extracted shape of W filter (used above Fl)
+%% X: extracted shape of X filter (used above Fl)
+%% Y: extracted shape of Y filter (used above Fl)
+%% Z: extracted shape of Z filter (used above Fl)
+%%    Z is optional, if not supplied X and Y are averaged to obtain Z
+%%
+%% 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/>.
+
+function [FIR] = sphear_design_fir_filters(F, n, maxboost, SR, W, X, Y, Z)
+  if ! exist("Z")
+    %% if Z was not supplied, just average X and Y data
+    Z = (X + Y) / 2;
+  end
+  %% make a 0 to 1 vector of relative frequencies
+  freqs = [0, (F/(SR/2))', 1]';
+  %% design the filters, we use twice the requested size as the
+  %% final filters will become 1/2 the size because of the minimal
+  %% phase calculation
+  IR(1,:) = fir2((n * 2) - 1, freqs, [W(1), W', W(end)]');
+  IR(2,:) = fir2((n * 2) - 1, freqs, [X(1), X', X(end)]');
+  IR(3,:) = fir2((n * 2) - 1, freqs, [Y(1), Y', Y(end)]');
+  IR(4,:) = fir2((n * 2) - 1, freqs, [Z(1), Z', Z(end)]');
+
+  for b = (1:4)
+    if (1)
+    %% try to get the filter to be close to minimum phase
+    %%
+    %% "Two alternative minimum-phase filters tested perceptually"
+    %% Robert Mores and Ralf Hendrych, AES Convention paper 9554
+    %%
+    %% number of samples in the sin transition
+    %%
+    %% 8 seems to be the sweet spot for a FIR filter size of 256,
+    %% 16 leads to bigger phase shift, 4 starts to impact negatively
+    %% the flattness of the frequency response
+    %%
+    M = size(IR)(2);
+    N = M/32;
+    R = 0.5 + (0.5 * sin(pi * ((1:(2*N)) - N) / (2 * N)));
+    ramp = cat(1, zeros(M/2 - N, 1), R', ones(M/2 - N, 1));
+    IRMP = IR(b,:)' .* ramp * 2;
+    %%
+    %% now we trim the filter and return one that is 1/2 the size requested
+    %% because it is now almost minimum phase and the first half is zeroed
+    %%
+    %% we also window the end of the impulse response... (unnecessary?)
+    %% T = 4;
+    %% tapper = cat(1, ones(M/2 - 2*T + 1, 1), (1.0 - (0.5 + (0.5 * sin(pi * ((1:(2*T)) - T) / (2 * T)))))');
+    %%
+    FIR(b,:) = IRMP(M/2 - N:M/2 - N + M/2 - 1);
+    else
+      size(IR(b,:))
+      FIR(b,:) = mps(fft(IR(b,:)));
+    end
+  end
+end
+

calibration/src/sphear_octa_af_calibrate.m

@@ -0,0 +1,83 @@
+%% -*- mode: octave -*-
+%%
+%% *SpHEAR Project
+%%
+%% sphear_af_octa_calibrate: read reference microphone signal, filter with
+%%                           convolution kernel created by DRC and filter
+%%                           all recorded signals, calibrated A format files
+%%                           are written to disk
+%%
+%% usage:
+%%   [] = sphear_af_octa_calibrate(measurements, keep, fade)
+%%
+%% keep: time to keep after main impulse in seconds
+%% fade: fade out time after keep in seconds
+%%
+%% 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/>.
+
+function [] = sphear_octa_af_calibrate(measurements, keep, fade);
+  %% 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
+  %%
+  %% read the convolution kernel
+  [ref, Fs] = wavread('../data/current/cal/drc-filter.wav');
+
+  %% window the DRC generated filter, eliminate the stuff before the main
+  %% peak and limit the length to the approximate useful length of the
+  %% original recorded impulse response
+  
+  [peaks, locs] = findpeaks(abs(ref));
+  [value, index] = max(peaks);
+  loc = locs(index);
+  
+  %% truncate before the main peak
+  b_len = 4;
+  before = tukeywin(2 * b_len, 0.9)(1:end/2) .* ref(loc-2*b_len:loc-b_len-1);
+  %% useful samples
+  samples = keep * Fs/2;
+  during = ref(loc-b_len:loc+samples-1);
+  %% fade out after useful samples
+  a_len = fade * Fs/2;
+  after = tukeywin(2 * a_len, 0.9)(end/2:end) .* ref(loc+samples:loc+samples+a_len);
+  %% and the whole thing is:
+  trunc = cat(1, before, during, after);
+
+  capnames = cellstr(['LF'; 'RF'; 'LB'; 'RB']);
+  octanames = cellstr(['C1'; 'C2'; 'C3'; 'C4'; 'C5'; 'C6'; 'C7'; 'C8']);
+  %% all requested measurements
+  for m = measurements
+    printf("processing measurement %d\n", m);
+    %% read in capsule signals and normalize it
+    for c = 1:4
+      %% first set of capsules (both sets have to have the same length)
+      sig = wavread(["../data/current/af/a-format1-" sprintf("%03d", m) "-" char(capnames(c)) ".wav"]);
+      wavwrite(sig, Fs, ["../data/current/af/a-format-" sprintf("%03d", 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"]);
+    end
+    for c = 1:4
+      %% second set of capsules
+      sig = wavread(["../data/current/af/a-format2-" sprintf("%03d", m) "-" char(capnames(c)) ".wav"]);
+      wavwrite(sig, Fs, ["../data/current/af/a-format-" sprintf("%03d", m) "-" char(octanames(c + 4)) ".wav"]);
+      cals = conv(trunc, sig);
+      cal_sig(:,c + 4) = conv(trunc, sig);
+      wavwrite(cal_sig(:,c + 4), Fs, ["../data/current/af/a-format-cal-" sprintf("%03d", m) "-" char(octanames(c + 4)) ".wav"]);
+    end
+    wavwrite(cal_sig, Fs, ["../data/current/af/a-format-cal-" sprintf("%03d", m) ".wav"]);
+  end
+end

calibration/src/sphear_plot_directions.m

@@ -0,0 +1,158 @@
+%% -*- mode: octave -*-
+%%
+%% *SpHEAR Project
+%%
+%% sphear_plot_directions: render amplitude response plots of B format signals
+%%                            version for current prototype
+%%
+%% usage:
+%%   [] = sphear_plot_directions(BF, fmin, fmax, ymin, ymax, Fs, Nfft, plot_on, smooth)
+%%
+%% fmin, fmax: range of frequencies to plot (Hz)
+%% ymin, ymax: vertical scale for plot (dB)
+%% Nfft: size of FFT
+%%
+%% optional parameters:
+%%
+%% plot_on: plot frequency response of components, default is "epsc"
+%%          can be one of 0, "eps", "epsc" (color), "png", "pngc" (color)
+%% smooth: number of samples to average, default is 0
+%%
+%% implied parameter:
+%%     files: "../data/current/bf-tetraproc/b-format-xxx.wav"
+%%     four channel B format, FUMA ordering and weights
+%%
+%% example:
+%%   sphear_plot_directions(lf_limit, hf_limit, -15, 3, render_fft_size, 1, 0);
+%%
+%% Copyright 2016, 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/>.
+
+function [] = sphear_plot_directions(BF, fmin, fmax, ymin, ymax, Fs, Nfft, plot_on, smooth)
+
+  %% no smoothing or plotting by default
+  if (!exist("smooth"))
+    smooth = 0;
+  end
+  if (!exist("plot_on"))
+    plot_on = "epsc";
+  else
+    if (plot_on == 1)
+      plot_on = "epsc";
+    end
+  end
+
+  %% plot if requested
+  if (plot_on != 0)
+    if (strcmp(plot_on, "epsc") || strcmp(plot_on, "pngc"))
+      %% color plot
+      wopts = 'k;W;';
+      xopts = 'r;X;';
+      yopts = 'b;Y;';
+      copt = '-color';
+      cext = "_color";
+    else
+      %% mono plot
+      wopts = '-';
+      xopts = '--';
+      yopts = ':';
+      copt = '-mono';
+      cext = "_mono";
+    end
+    f = figure;
+    hold on;
+    grid on;
+    axis([fmin, fmax, ymin, ymax]);
+    tics('y', [6, 3, 0, -3, -6, -9, -12, -15, -20, -25, -30]);
+    tics('x', [300, 500, 1000, 1500, 2000, 3000, 5000, 7000, 10000, 15000, 20000]);
+    ca=get (gcf, "currentaxes");
+    set(ca,"fontweight","bold","linewidth",1.5);
+
+    %% measurements, first four in cardinal directions
+    measurements = [1, 5, 9, 13, 3, 7, 11, 15];
+
+    %% which components to plot, we ommit the ones that are nulled
+    %% 0: do not plot, 1: cardinal, 2: diagonal plot
+    m_W = [1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0];
+    m_X = [1, 0, 2, 0, 0, 0, 2, 0, 1, 0, 2, 0, 0, 0, 2, 0];
+    m_Y = [0, 0, 2, 0, 1, 0, 2, 0, 0, 0, 2, 0, 1, 0, 2, 0];
+    m_Z = [1, 0, 2, 0, 1, 0, 2, 0, 1, 0, 2, 0, 1, 0, 2, 0];
+    
+    idx = 1;
+    for m = measurements
+      %% calculate B format frequency responses
+      for b = 1:4
+	if (smooth > 0)
+	  bfft(idx,b,:) = fastsmooth(10*log10(abs(fft(BF(m, b, :), Nfft))), smooth, 2, 1)(1:end/2);
+	else
+	  bfft(idx,b,:) = 10*log10(abs(fft(BF(m, b, :), Nfft)))(1:end/2);
+	end
+	idx = idx + 1;
+      end
+    end
+
+    %% scale everything to average of W (1st measurement)
+    fft_avg = mean(bfft(1,1,round(end/4):round(end/2)));
+
+    %% now do the plots
+    idx = 1;
+    W_offset = +0;
+    X_offset = -6;
+    Y_offset = -6;
+    fr = (1:Nfft/2) * Fs/Nfft;
+    for m = measurements
+      for b = 1:4
+	if (m_W(m) > 0 && b == 1)
+	  if (m_W(m) > 1)
+	    pW = semilogx(fr, bfft(idx,b,:) - fft_avg + W_offset, wopts, "linewidth", 1.5, "color", 'k', "linestyle", "--");
+	  else
+	    pW = semilogx(fr, bfft(idx,b,:) - fft_avg + W_offset, wopts, "linewidth", 1.5, "color", 'k');
+	  end
+	end
+	if (m_X(m) > 0 && b == 2)
+	  if (m_X(m) > 1)
+	    pXd = semilogx(fr, bfft(idx,b,:) - fft_avg + X_offset, xopts, "linewidth", 1.5, "color", "r", "linestyle", "--");
+	  else
+	    pXp = semilogx(fr, bfft(idx,b,:) - fft_avg + X_offset, xopts, "linewidth", 1.5, "color", 'r');
+	  end
+	end
+	if (m_Y(m) > 0 && b == 3)
+	  if (m_Y(m) > 1)
+	    pYd = semilogx(fr, bfft(idx,b,:) - fft_avg + Y_offset, yopts, "linewidth", 1.5, "color", "b", "linestyle", "--");
+	  else
+	    pYp = semilogx(fr, bfft(idx,b,:) - fft_avg + Y_offset, yopts, "linewidth", 1.5, "color", 'b');
+	  end
+	end
+	idx = idx + 1;
+      end
+    end
+    l = legend([pW, pXp, pXd, pYp, pYd],
+	       {sprintf("W principal directions (%.1fdB offset)", W_offset),
+		sprintf("X principal directions (%.1fdB offset)", X_offset),
+		sprintf("X diagonal directions (%.1fdB offset)", Y_offset),
+		sprintf("Y principal directions (%.1fdB offset)", X_offset),
+		sprintf("Y diagonal directions (%.1fdB offset)", Y_offset),
+	       });
+    legend(l, 'location', 'southwest');
+    if (strcmp(plot_on, 'eps') || strcmp(plot_on, 'epsc') || plot_on == 1)
+      print(f, '-deps', '-tight', '-S1920,1080', copt, ["../figures/b_format_sphear3.eps"]);
+    else
+      print(f, '-deps', copt, ["../figures/b_format_sphear3.eps"]);
+    end
+
+  end
+end
+

calibration/src/sphear_write_simple_a2b_encoder.m

@@ -0,0 +1,94 @@
+%% -*- mode: octave -*-
+%%
+%% *SpHEAR Project
+%%
+%% sphear_write_a2b_encoder: write simple Faust code for an A format to B format converter
+%%                           (static A to B matrix plus 4 FIR filters, one of each B component)
+%%
+%% usage:
+%%   [] = sphear_write_simple_a2b_encoder(A2B, FIR)
+%%
+%% A2B: static A to B matrix
+%% FIR: FIR filter matrix
+%%
+%% Copyright 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/>.
+
+function [] = sphear_write_simple_a2b_encoder(A2B, FIR)
+
+  filename = '../data/current/cal/simple_a2b_encoder.dsp';
+  [file, msg] = fopen(filename,'w');
+  if msg
+    error([msg ': ' filename]);
+  end
+
+  %% write beginning of file...
+  fprintf(file, 'import("stdfaust.lib");\n');
+  fprintf(file, '\n');
+
+  bnames = cellstr(['w'; 'x'; 'y'; 'z']);
+
+  %% write static A to B matrix rows
+  fuma = [sqrt(2)/2, 1, 1, 1];
+  for a = 1:size(A2B,1)
+    for b = 1:4
+      fprintf(file, ["a" sprintf("%d", a) char(bnames(b)) " = "]);
+      write_faust_vector(file, fuma(b) * A2B(a, b));
+      fprintf(file, ';\n\n');
+    end
+  end
+  
+  %% write four FIR filter coefficients
+  for b = 1:4
+    bf_filter = reshape(FIR(b, :), 1, []);
+    fprintf(file, [char(bnames(b)) "_fir = "]);
+    write_faust_vector(file, bf_filter);
+    fprintf(file, ';\n\n');
+  end
+  
+  %% be done with it...
+  if (size(A2B,1) == 4)
+    fprintf(file, 'process = _,_,_,_ <: w*wgain, x*gain, y*gain, z*gain \n');
+  else
+    %% assume 8 input channels for now
+    fprintf(file, 'process = _,_,_,_,_,_,_,_ <: w*wgain, x*gain, y*gain, z*gain \n');
+  end
+  fprintf(file, 'with {\n');
+  if (size(A2B,1) == 4)
+    fprintf(file, '     w = *(a1w), *(a2w), *(a3w), *(a4w) :> fi.fir(w_fir) ;\n');
+    fprintf(file, '     x = *(a1x), *(a2x), *(a3x), *(a4x) :> fi.fir(x_fir) ;\n');
+    fprintf(file, '     y = *(a1y), *(a2y), *(a3y), *(a4y) :> fi.fir(y_fir) ;\n');
+    fprintf(file, '     z = *(a1z), *(a2z), *(a3z), *(a4z) :> fi.fir(z_fir) ;\n');
+  else
+    fprintf(file, '     w = *(a1w), *(a2w), *(a3w), *(a4w), *(a5w), *(a6w), *(a7w), *(a8w) :> fi.fir(w_fir) ;\n');
+    fprintf(file, '     x = *(a1x), *(a2x), *(a3x), *(a4x), *(a5x), *(a6x), *(a7x), *(a8x) :> fi.fir(x_fir) ;\n');
+    fprintf(file, '     y = *(a1y), *(a2y), *(a3y), *(a4y), *(a5y), *(a6y), *(a7y), *(a8y) :> fi.fir(y_fir) ;\n');
+    fprintf(file, '     z = *(a1z), *(a2z), *(a3z), *(a4z), *(a5z), *(a6z), *(a7z), *(a8z) :> fi.fir(z_fir) ;\n');
+  end
+  
+  fprintf(file, '     gain = hslider("XYZ",1,0,1,0.01);\n');
+  fprintf(file, '     wgain = hslider("W",1,0,1,0.01);\n');
+  fprintf(file, '};\n');
+  fprintf(file, '\n');
+  fclose(file);
+end
+
+function [] = write_faust_vector(file, v)
+  v = v(:);
+  fprintf(file, '(');
+  fprintf(file, '% 14.10g,', v(1:end-1));
+  fprintf(file, '% 14.10g)', v(end));
+end