add sphear_fir_render_directions, trying to get a graph similar to Richard Lee's analysis

sphear_fir_render_directions.m

@@ -0,0 +1,109 @@
+%%
+%% *SpHEAR Project
+%% sphear_fir_render_directions: render impulse responses throught the FIR filter A2B matrix
+%%
+%% usage:
+%%   [] = sphear_fir_render_directions(FIR, fmin, fmax, N, Npre, Fs, 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
+%% Fs: sampling rate
+%%
+%% 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 [BFF] = sphear_fir_render_directions(FIR, fmin, fmax, N, Npre, Fs, smooth)
+
+  [ref, Fs] = wavread('../cal/2016.04.20/export/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);
+  
+  capnames = cellstr(['LF'; 'RF'; 'LB'; 'RB']);
+
+  minindex = N;
+  maxpeak = 0;
+  
+  Nfft = 4096;
+
+  f = figure;
+  hold on;
+  grid on;
+  axis([fmin, fmax, -30, 0]);
+
+  %% title(sprintf("azimuth = %05.1f degrees", angle));
+  tics('y', [6, 0, -3, -6, -9, -12, -15, -20, -25, -30]);
+  tics('x', [300, 500, 1000, 2000, 3000, 5000, 7000, 10000, 15000, 20000]);
+
+  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);
+
+    %% calculate B format frequency responses
+    for b = 1:4
+      bfft(idx,b,:) = fastsmooth(20*log10(abs(fft(bf(idx,b,:), Nfft))), smooth, 2, 1);
+    end
+    fr = (1:Nfft/2) * Fs/Nfft;
+    semilogx(fr, bfft(idx,1,1:end/2), 'k');
+    if (idx > 4)
+      col = 'b';
+    else
+      col = 'r';
+    end
+    if (rem(idx, 2) == 1)
+      semilogx(fr, bfft(idx,2,1:end/2), col);
+    else
+      semilogx(fr, bfft(idx,3,1:end/2), col);
+    end
+
+    idx = idx + 1;
+  end
+
+  print(f, '-deps','-color',["figures/sphear_directions.eps"]);
+end
+
+
+function [bf] = sphear_a_to_b(m, FIR, range, N);
+  capnames = cellstr(['LF'; 'RF'; 'LB'; 'RB']);
+  %% read in capsule signals and apply FIR filter
+  for c = 1:4
+    cs(c,:) = wavread(["../cal/2016.04.20/export/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,:), N);
+    end
+  end
+  %% add filtered capsule output to b format components
+  for b = 1:4
+    bf(b,:) = bs(1,b,:) + bs(2,b,:) + bs(3,b,:) + bs(4,b,:);
+  end
+end
+
+
+function [i, v] = findmaxpeak(signal);
+  [peaks, locs] = findpeaks(signal, "DoubleSided");
+  [~, index] = max(peaks);
+  i = locs(index);
+  v = peaks(index);
+end