add functions to create polar plots for encoder created by Fons

record_fons_b_format

@@ -0,0 +1,18 @@
+#!/bin/bash
+#
+
+for m in `seq 1 16`; do
+    n=`printf "%03d" $m`
+    echo "==== measurement $m ($n)..."
+    # start background capture
+    ecasound -t:6 \
+	     -i "jack_multi,tetraproc:B-form.W,tetraproc:B-form.X,tetraproc:B-form.Y,tetraproc:B-form.Z" \
+	     -f:f32,4,48000 -o:sndfile,"../cal/2016.04.20/export/b-format-fons-rendered-$n.wav" &
+    # play back A format signal
+    ecasound -f:s16_le,4,48000 -i:sndfile,../cal/2016.04.20/export/a-format-fons-$n.wav \
+	     -o "jack_multi,tetraproc:A-form.1(LF),tetraproc:A-form.2(RF),tetraproc:A-form.3(LB),tetraproc:A-form.4(RB)" \
+	     -eadb:0
+    # wait for it...
+    wait
+done
+

sphear_a2b_fons_all.m

@@ -0,0 +1,62 @@
+%%
+%% *SpHEAR Project
+%%
+%% sphear_a2b_all: calculate A2B matrixes for frequency ranges
+%%
+%% usage:
+%%   [F FE A2BF] = sphear_a2b_all(m, fmin, fmax, Noct, Fs, Nsamp, Npre, 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
+%% Fs: sampling rate
+%% Nsamp: time domain samples to use
+%% Npre: samples to use before the main impulse peak
+%%
+%% 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 [B] = sphear_a2b_fons_all(fmin, fmax, Noct, Nfft, Fs, plt)
+
+  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_bf_power_fons(600, 1200, Nfft, Fs, 1);
+  norm = mean(ref(:,1));
+
+  for i = (1:size(fc))
+    B(:,:,i) = sphear_bf_power_fons(fe(i), fe(i + 1), Nfft, Fs, plotz) / norm;
+    if (plt == 1)
+      sphear_bf_plot(B(:,:,i), fe(i), fe(i + 1), 3.5);
+    end
+  end
+
+end
+

sphear_bf_plot_fons.m

@@ -0,0 +1,92 @@
+%%
+%% *SpHEAR Project
+%%
+%% sphear_bf_plot: plot the polar response of a B format signal
+%%
+%% usage:
+%%   [BF] = sphear_bf_plot(BF, fmin, fmax)
+%%
+%% BF: B format signal array
+%% fmin, fmax: range of frequencies to consider (to title plot)
+%%
+%% Copyright Fernando Lopez-Lezcano 2016
+%% 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 [BF] = sphear_bf_plot_fons(BF, fmin, fmax, rmax)
+
+  if (!exist("plot_on"))
+    plot_on = "epsc";
+  else
+    if (plot_on == 1)
+      plot_on = "epsc"
+    end
+  end
+
+  %% to close the polar plots we repeat the first point
+  m = size(BF, 1);
+  n = 32;
+  BF(m + 1,:) = BF(1,:);
+  theta = (0:m)*(2*pi/m);
+  thetaf = (0:n)*(2*pi/n);
+  
+  %% plot it as a function of angle
+  bname = cellstr(['W'; 'X'; 'Y'; 'Z']);
+
+  f = figure;
+  for b = 1:4
+    subplot(2, 2, b);
+    hold on;
+    grid off;
+    axis square;
+    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"]);
+
+    if (strcmp(plot_on, "epsc") || strcmp(plot_on, "pngc"))
+      %% color plot
+      dopts = '.-r';
+      ropts = '-b';
+      copt = '-color';
+      cext = "_color";
+    else
+      %% mono plot
+      dopts = '.-';
+      ropts = '-';
+      copt = '-mono';
+      cext = "_mono";
+    end
+
+    grey = [0.1, 0.1, 0.1];
+    %% radius circles
+    for rad = 0.25:0.25:rmax
+      [xr, yr] = pol2cart(thetaf', ones(n + 1,1) * rad);
+      plot(xr, yr, ':', 'Color', grey);
+    end
+    %% radial lines
+    for ang = theta
+      plot([0 rmax*cos(ang)], [0 rmax*sin(ang)], ':', 'Color', grey);
+    end
+
+    [xr, yr] = pol2cart(thetaf', ones(n + 1,1));
+    plot(xr, yr, ropts);
+    [x, y] = pol2cart(theta', abs(BF(:,b)));
+    plot(x, y, dopts, 'linewidth', 1.5, 'markersize', 6);
+    
+  end
+  print(f, '-deps','-tight',copt,["figures/sphear_polar_multi_" sprintf("%05d", fmin) "_" sprintf("%05d", fmax) cext ".eps"]);
+end
+

sphear_bf_power_fons.m

@@ -0,0 +1,63 @@
+%%
+%% *SpHEAR Project
+%%
+%% sphear_a_power: return array of power for all capsules in a range of frequencies
+%%
+%% usage:
+%%   [R] = sphear_a_power(fmin, fmax, N, Npre)
+%%
+%% fmin, fmax: range of frequencies to consider
+%% Npre: number of samples before the main impulse
+%% N: total number of samples to analyze
+%%
+%% The names of reference microphone impulse response and capsule measurements are
+%% currently hardwired into this file, this needs to be fixed and turned into
+%% input parameters
+%%
+%% Copyright Fernando Lopez-Lezcano 2016
+%% 
+
+function [B] = sphear_bf_power_fons(fmin, fmax, N, Fs, plotz)
+  pkg load signal;
+
+  %% indexes in array for frequency bins
+  imin = round((N/2) * fmin / (Fs/2));
+  imax = round((N/2) * fmax / (Fs/2));
+
+  r_start = 0;
+  r_end = 0;
+  
+  %% all measurements
+  for m = 1:16
+    %% all b format components, 4 channel soundfile
+    bsignal = wavread(["../cal/2016.04.20/export/b-format-fons-rendered-" sprintf("%03.0f", m) "-16.wav"]);
+
+    %% transform from FUMA to SN3D (the encoder is outputing FUMA)
+    %% bsignal(:,1) = bsignal(:,1) * (sqrt(2)/2);
+	
+    %% find first peak, read samples around it
+    Npre = 64;
+    Nbuf = 256;
+
+    [index, value] = findmaxpeak(bsignal(:,1));
+    r_start = index - Npre;
+    r_end = r_start + Nbuf -1;
+    range = r_start:r_end;
+    printf("%d: samples = %d -> %d (%d), imin:imax %d:%d\n", m, r_start, r_end, r_end-r_start, imin, imax);
+
+    bfft = fft(bsignal(range,:), N);
+    if (plotz == 0)
+      %% calculate a forward facing cardiod
+      cfft = fft(bsignal(range,1) + 1.0 * bsignal(range,2), N);
+    end
+    %% return total power of signal in frequency range
+    scaler = 8;
+    B(m, :) = sqrt(sum(bfft(imin:imax,:).*conj(bfft(imin:imax,:)))/size(bfft(imin:imax,:),1)) / scaler;
+    %% make a cardiod if we don't want Z
+    if (plotz == 0)
+      B(m, 4) = sqrt(sum(cfft(imin:imax).*conj(cfft(imin:imax)))/size(cfft(imin:imax),1)) / (2 * scaler);
+    end
+  end
+  
+end
+