calibration: add sphear_bf_plot_directions: plot principal and diagonal directions for signals recorded from tetraproc and using the calculated FIR filter matrix; record_tetraproc_b_format: pyecasound script to play calibrated IRs through tetraproc and record the B format output to files

calibration/src/record_tetraproc_b_format

@@ -0,0 +1,109 @@
+#!/bin/python
+#
+# PyEcasound Control Script
+#
+# Play back A format files (already trimmed and windowed) from disk to
+# the inputs of TetraProc, and at the same time record B format data from
+# the output of TetraProc to a disk file.
+#
+# As ecasound does not recognize the input soundfiles written by octave,
+# they are converted first to 24 bit PCM using sndfile-convert
+#
+#
+# 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/>.
+
+import time
+from ecacontrol import *
+from subprocess import *
+
+e = ECA_CONTROL_INTERFACE()
+
+e.command("cs-add tetraproc_process")
+
+status = e.command("c-add AF_play")
+status = e.command("c-add BF_record")
+print(status)
+
+## new chain for A-format playback
+
+print(e.command("c-select AF_play"))
+print("== starting jack output...")
+print(e.command("cs-set-audio-format ,4,48000"))
+print(e.command("ao-add jack_multi,tetraproc:A-form.1(LF),tetraproc:A-form.2(RF),tetraproc:A-form.3(LB),tetraproc:A-form.4(RB)"))
+
+## new chain for B-format recording
+
+print(e.command("c-select BF_record"))
+print("== starting jack input...")
+print(e.command("cs-set-audio-format ,4,48000"))
+print(e.command("ai-add jack_multi,tetraproc:B-form.W,tetraproc:B-form.X,tetraproc:B-form.Y,tetraproc:B-form.Z"))
+
+## start processing
+
+for measurement in range(1, 17):
+    m = "%0.3d" % (measurement)
+    print("==== processing measurement %d" % (measurement))
+
+    ## convert soundfile to format needed by ecasound
+    os.system("sndfile-convert -pcm24 ../data/current/af/a-format-cal-" + m + ".wav " +
+              " ../data/current/af/a-format-cal-" + m + "-pcm24.wav ")
+    
+    ## set up file playbak
+    print(e.command("c-select AF_play"))
+    print("== reading soundfile...")
+    print(e.command("ai-add ../data/current/af/a-format-cal-" + m + "-pcm24.wav"))
+
+    ## set up file recording
+    print(e.command("c-select BF_record"))
+    print("== writing soundfile...")
+    print(e.command("ao-add ../data/current/bf-tetraproc/b-format-" + m + ".wav"))
+
+    ## start jack...
+    print(e.command("engine-launch"))
+    time.sleep(0.25)
+
+    print(e.command("start"))
+    print("== waiting...")
+    time.sleep(1.0)
+    print("== woke up")
+
+    e.command("stop")
+    print("== stopped")
+
+    ## remove old chains
+    # (this prints an error but seems to work anyway)
+    print("== remove objects")
+    print(e.command("ao-select ../data/current/bf-tetraproc/b-format-" + m + ".wav"))
+    print(e.command("ao-remove"))
+
+    ## wait a little bit...
+    time.sleep(0.5)
+
+    #print(e.command("aio-status"))
+    
+## done with processing
+
+e.command("cs-disconnect")
+print("== disconnected!")
+
+## remove all temporary files
+for measurement in range(1, 17):
+    m = "%0.3d" % (measurement)
+    print("==== removing %d" % (measurement))
+    os.system("rm -f ../data/current/af/a-format-" + m + "-pcm24.wav")
+
+print("== done")

calibration/src/sphear_bf_plot_directions.m

@@ -0,0 +1,165 @@
+%% -*- mode: octave -*-
+%%
+%% *SpHEAR Project
+%%
+%% sphear_bf_plot_directions: render amplitude response plots of B format signals
+%%                            version for current prototype
+%%
+%% usage:
+%%   [] = sphear_bf_plot_directions(fmin, fmax, ymin, ymax, 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_bf_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_bf_plot_directions(fmin, fmax, ymin, ymax, 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
+      %% read rendered four channel B format measurement
+      %% trim to just the needed part
+      range = [512:4096];
+      [bsig, Fs] = wavread(["../data/current/bf-tetraproc/b-format-" sprintf("%03d", m(1)) ".wav"]);
+      %% convert tetraproc FUMA output to SN3D
+      bsig(:,1) = bsig(:,1) / 0.707;
+      %% calculate B format frequency responses
+      for b = 1:4
+	bsig_fft = fft(bsig(range,b), Nfft);
+	if (smooth > 0)
+	  bfft(idx,b,:) = fastsmooth(10*log10(abs(bsig_fft)), smooth, 2, 1)(1:end/2);
+	else
+	  bfft(idx,b,:) = 10*log10(abs(bsig_fft))(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
+