cal: add encoder design example for the Octathingy (work in progress)

calibration/src/sphear_octa_design_example.m

@@ -0,0 +1,305 @@
+%% -*- mode: octave -*-
+%%
+%% *SpHEAR Project
+%%
+%% Design a second order A2B encoder for the Octathingy
+%%
+%% Eight capsule equalization filters, static A to B matrix
+%% and eight filters for the B-format components
+%%
+%% Copyright 2016-2018, 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/>.
+
+%% root of the data files, assume we are running in the src/ directory
+DATA = "../data/";
+
+%% default for location of current measurements
+CURRENT = "current/";
+
+%% data is assumed to be in ../data/current/ which is usually a symbolic
+%% link to the real measurement directory
+
+%% before we start the design process we need to equalize the measurements
+%% to take into account the frequency response of the speaker used. The first
+%% step is to use the reference microphone measurement to design a filter that
+%% inverts the frequency response of the speaker. We use DRC to do this. 
+
+%% calculate speaker calibration filter from reference microphone
+%% (shell code, uses DRC for filter calculation, edit script to change
+%% the DRC preset to use, the default uses the "minimal" preset)
+
+%% ./sphear_make_reference_filter emm-6-5850
+
+%% at this point we should have the raw impulse responses and a calibration
+%% filter in ../data/current/cal/drc-filter.wav (see the script for details
+%% on how this is done). We process all impulse responses through the filter
+%% and write calibrated impulse responses that are used for the rest of the
+%% design process. 
+
+%% read raw impulse response data and generate calibrated impulse responses
+%% (this is a version for the Octathingy which relies on recording two four
+%% channel impulse responses in Aliki - keep in mind that both recordings
+%% have to be processed and trimmed alike, otherwise there will be a mismatch
+%% in arrival times of the impulses)
+
+sphear_octa_af_calibrate((1:16), 0.0045, 0.001);
+sphear_octa_af_calibrate((1:16), 0.0045, 0.001, "up-");
+sphear_octa_af_calibrate((1:16), 0.0045, 0.001, "down-");
+
+%% read in calibrated signal (reads all three layers)
+
+[AF, FS] = sphear_af_read((1:16));
+
+%% capsule frequency response as a function of incident angle
+
+sphear_capsule_plot_frequency_response(AF, [1:16], 300, 20000, -60, -10, 48000, 8192);
+
+%% same for the other two layers
+%% up
+sphear_capsule_plot_frequency_response(AF, [17:32], 300, 20000, -60, -10, 48000, 8192);
+%% down: we see the rig interference in the low frequency ripple
+%%       we should probably trim it to remove that part of the IR...
+sphear_capsule_plot_frequency_response(AF, [33:48], 300, 20000, -60, -10, 48000, 8192);
+
+%% frequency response of all capsules for each incident angle
+
+sphear_af_plot_frequency_response(AF, 300, 20000, -60, -10, 48000, 8192);
+
+%% calculate signal power in reference frequency interval for all
+%% measurements and capsules
+
+[R] = sphear_signal_power(AF, FS, 800, 2000, 8192);
+
+%% check that capsule signals are reasonable
+
+sphear_plot_r(R(1:16,:), 0);
+sphear_plot_r(R(17:32,:), 0);
+sphear_plot_r(R(33:48,:), 0);
+
+%% as the capsules are tilted with respect to the horizontal plane we never
+%% hit the null in the back, the effect is more pronounced for the "up" and
+%% "down" measurements
+
+%% now we need to delay the signals for all capsules for each measurement
+%% so that all individual measurements are aligned in time.
+%%
+%% we start by calculating the average delay for all capsules in each
+%% measurement (delays are relative and calculate taking the first capsule
+%% of the first measurement as a reference)
+
+%% using cross correlation (deprecated)
+%% [DELAYS] = sphear_signal_delta_arrival_times(AF);
+%% using parabolic interpolation (use this one)
+[DELAYS DELTA ALL] = sphear_signal_arrival_times(AF);
+
+%% and then we delay the signals with a fractional sync delay
+
+[AFC] = sphear_fix_wobble(AF, DELAYS);
+
+%% this will fix the wobble between measurements as the array is not
+%% perfectly centered when rotating in the measuring rig, it will
+%% also fix any inconsistencies between different measurement layers
+%% if the rig has to be moved and re-aligned between layers
+
+%% another run tells us if the signals are now aligned...
+
+[ADELAYS ADELTA AALL] = sphear_signal_arrival_times(AFC);
+
+%% and a plot to see how delays were fixed
+
+plot(DELTA, "b.-");
+hold on
+grid on
+plot(ADELTA, "r.-");
+plot([16,16],[0,max(DELTA)], "k");
+plot([32,32],[0,max(DELTA)], "k");
+hold off
+
+%% the result is not perfect... why? Most likely because the peaks
+%% are not always chosen correctly, see below for one example
+
+%% the alignment depends on finding the arrival time of the impulse for each
+%% measurement, the following code plots that so that we can see if
+%% things are happening the way they should
+
+%% first measurement, first capsule
+plot(abs(AF(1,1,1:128)), "b.-");
+hold on
+grid on
+%% find peaks
+[amps, locs, widths, residue] = sphear_findpeaks(abs(reshape(AF(1,1,1:128), 1, [])), 10);
+%% add a line for the peak
+plot([locs(1), locs(1)], [0.0, amps(1)], "k.-");
+%% first measurement, fifth capsule (back)
+plot(abs(AF(1,5,1:128)), "r.-");
+%% find peaks
+[amps, locs, widths, residue] = sphear_findpeaks(abs(reshape(AF(1,5,1:128), 1, [])), 10);
+%% add a line for the peak
+plot([locs(1), locs(1)], [0.0, amps(1)], "k.-");
+%% plot the bipolar signals as well
+plot(AF(1,1,1:128), "b.--");
+plot(AF(1,5,1:128), "r.--");
+hold off
+
+%% plot a single measurement (all capsules) [needs to be fixed]]
+%% THETA=(2*pi/16:-2*pi/8:-2*pi+2*pi/8)*180/pi
+%% polar(THETA,ALL(2,:) - mean(ALL(2,:)), "b.-");
+
+%% calculate full filter set
+
+[R] = sphear_signal_power(AFC, FS, 800, 2000, 8192);
+
+sphear_plot_r(R(1:16,:), 0);
+sphear_plot_r(R(17:32,:), 0);
+sphear_plot_r(R(33:48,:), 0);
+
+%% create X matrix for three layers
+
+M = sphear_X_matrix(16, 43, -43);
+
+%% X matrix without S/T columns
+
+MH = M([1:16],:);
+
+RH = R([1:16],:);
+
+%% horizontal only A2B for WXYUV
+%% we use the built-in X matrix in A2B_matrix.m
+
+[A2BH, PV, COND] = sphear_calculate_a2b(RH);
+
+%% convert A to B
+
+[BFH] = sphear_af_to_bf(AFC([1:16],:,:), A2BH);
+
+[BH, FC, FE] = sphear_signal_power_all(BFH, FS, 300, 20000, 800, 2000, 16, 8192);
+
+[W X Y Z S T U V] = sphear_bf_extract_filter_shapes(BH, 0);
+
+sphear_bf_plot_filters(FC, W, X, Y, Z, S, T, U, V);
+
+[FIRH IR MPS] = sphear_design_fir_filters(FC, 512, 48000, W, X, Y, Z, S, T, U, V, 400, 18, -10);
+
+[BFH, BFFH] = sphear_af_to_bf(AFC([1:16],:,:), A2BH, FIRH, 8192);
+
+[A2BH XY_gain ST_gain UV_gain] = sphear_bf_normalize_gain(A2BH, BFFH, FS, 800, 2000, 8192);
+
+[BFH, BFFH] = sphear_af_to_bf(AFC([1:16],:,:), A2BH, FIRH, 8192);
+
+sphear_bf_plot_frequency_response(BFFH, 600, 1200, 300, 20000, -30, 6, 48000, 8192);
+
+%% some polar plots
+
+sphear_bf_plot_polar_response(BFFH, 600, 1200, 300, 20000, 3, 8192, FS, 0, 0, 4);
+sphear_bf_plot_polar_response(BFFH, 600, 1200, 300, 20000, 3, 8192, FS, 0, 0, 12);
+sphear_bf_plot_polar_response(BFFH, 600, 1200, 300, 20000, 3, 8192, FS, 0, 0, 14);
+sphear_bf_plot_polar_response(BFFH, 600, 1200, 300, 20000, 3, 8192, FS, 0, 0, 18);
+
+%% bad performance at the top end of the spectrum for the 1st order signals... (15KHz and up)
+
+%% now calculate with the full 3d data set
+
+[A2B, PV, COND] = sphear_calculate_a2b(R, M);
+
+[BF] = sphear_af_to_bf(AFC, A2B);
+
+[B, FC, FE] = sphear_signal_power_all(BF, FS, 300, 20000, 800, 2000, 16, 8192);
+
+[W X Y Z S T U V] = sphear_bf_extract_filter_shapes(B, 0);
+
+[FIR IR MPS] = sphear_design_fir_filters(FC, 512, 48000, W, X, Y, Z, S, T, U, V, 400, 18, -10);
+
+[FIR IR MPS] = sphear_design_fir_filters(FC, 512, 48000, W, X, Y, Z, S, T, U, V, 0, 18, -50);
+
+[BF, BFF] = sphear_af_to_bf(AFC, A2B, FIR, 8192);
+
+[A2B XY_gain S_gain T_gain U_gain V_gain] = sphear_bf_normalize_gain(A2B, BFF, FS, 800, 2000, 8192);
+
+[BF, BFF] = sphear_af_to_bf(AFC, A2B, FIR, 8192);
+
+sphear_bf_plot_filters(FC, W, X, Y, Z, S, T, U, V);
+
+sphear_bf_plot_frequency_response(BFF(1:16,:,:), 600, 1200, 300, 20000, -30, 6, 48000, 8192);
+sphear_bf_plot_frequency_response(BFF(17:32,:,:), 600, 1200, 300, 20000, -30, 6, 48000, 8192);
+sphear_bf_plot_frequency_response(BFF(33:48,:,:), 600, 1200, 300, 20000, -30, 6, 48000, 8192);
+
+%% horizontal plane performance is pretty good, but up and down are bad at high frequencies
+%% second order components (specially S/T) not so good
+
+%%%%%%%%%
+
+%% try to equalize the capsules first
+
+[AH, FC, FE] = sphear_signal_power_all(AFC, FS, 300, 20000, 800, 2000, 16, 8192);
+
+%% how do we derive filters?
+%%
+%% we should try to equalize for frontal incident measurements, for example
+%% average the on axis plus the two measurements to the sides of the on-axis
+%% (in the horizontal plane, which catches the capsule slightly "sideways")
+
+[FSHAPE] = sphear_extract_capsule_filter_shapes(AH);
+
+%% design the filters
+
+[FIR IR MPS] = sphear_design_capsule_fir_filters(FC, 512, 48000, FSHAPE);
+
+%% filter the signals so they are now "flat"
+
+[AFF] = sphear_filter_capsule_signals(AFC, FIR, 8192);
+
+%% only the two adjacent measurements are no good, averaging all of them is also
+%% not the best - 
+
+sphear_capsule_plot_frequency_response(AFF, [1:16], 300, 20000, -60, -10, 48000, 8192);
+
+sphear_capsule_plot_frequency_response(AFC, [1:16], 300, 20000, -60, -10, 48000, 8192);
+sphear_capsule_plot_frequency_response(AFC, [17:32], 300, 20000, -60, -10, 48000, 8192);
+sphear_capsule_plot_frequency_response(AFC, [33:48], 300, 20000, -60, -10, 48000, 8192);
+
+%% what happens if we design the matrix...
+
+[R] = sphear_signal_power(AFF, FS, 800, 2000, 8192);
+
+[A2B, PV, COND] = sphear_calculate_a2b(R, M);
+
+[BF] = sphear_af_to_bf(AFF, A2B);
+
+[B, FC, FE] = sphear_signal_power_all(BF, FS, 300, 20000, 800, 2000, 16, 8192);
+
+[W X Y Z S T U V] = sphear_bf_extract_filter_shapes(B, 0);
+
+[FIR IR MPS] = sphear_design_fir_filters(FC, 512, 48000, W, X, Y, Z, S, T, U, V, 400, 18, -10);
+
+[FIR IR MPS] = sphear_design_fir_filters(FC, 512, 48000, W, X, Y, Z, S, T, U, V, 0, 18, -50);
+
+[FIR IR MPS] = sphear_design_fir_filters(FC, 512, 48000, W, X, Y, Z, S, T, U, V, 0);
+
+[BF, BFF] = sphear_af_to_bf(AFF, A2B, FIR, 8192);
+
+[A2B XY_gain S_gain T_gain U_gain V_gain] = sphear_bf_normalize_gain(A2B, BFF, FS, 800, 2000, 8192);
+
+[BF, BFF] = sphear_af_to_bf(AFF, A2B, FIR, 8192);
+
+sphear_bf_plot_filters(FC, W, X, Y, Z, S, T, U, V);
+
+sphear_bf_plot_frequency_response(BFF(1:16,:,:), 600, 1200, 300, 20000, -30, 6, 48000, 8192);
+
+
+sphear_bf_plot_polar_response(BFF(1:16,:,:), 600, 1200, 300, 20000, 3, 8192, FS, 0, 0, 4);
+sphear_bf_plot_polar_response(BFF, 600, 1200, 300, 20000, 3, 8192, FS, 0, 0, 12);
+sphear_bf_plot_polar_response(BFF, 600, 1200, 300, 20000, 3, 8192, FS, 0, 0, 14);
+sphear_bf_plot_polar_response(BFF, 600, 1200, 300, 20000, 3, 8192, FS, 0, 0, 18);