add sphear_fir_render.m, render impulse responses throught the FIR filter A2B matrix

sphear_fir_render.m

@@ -0,0 +1,107 @@
+%%
+%% *SpHEAR Project
+%% sphear_fir_render: render impulse responses throught the FIR filter A2B matrix
+%%
+%% usage:
+%%   [] = sphear_fir_render(FIR, fmin, fmax, N, Npre)
+%%
+%% 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
+%%
+%% 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_fir_render(FIR, fmin, fmax, N, Npre)
+
+  [ref, Fs] = wavread('../cal/2016.04.20/export/reference-001-1.wav');
+  %% find highest peak location and define sample range
+  [peaks, locs] = findpeaks(abs(ref));
+  [~, index] = max(peaks);
+  start = locs(index) - Npre;
+  range = start:(start + N - 1);
+  ref_fft = fft(ref(range), N);
+  
+  capnames = cellstr(['LF'; 'RF'; 'LB'; 'RB']);
+
+  raw = 0;
+
+  for m = 1:16
+    c1 = wavread(["../cal/2016.04.20/export/a-format-0" sprintf("%02.0f", m) "-" char(capnames(1)) ".wav"])(range);
+    w1 = fftfilt(reshape(FIR(1,1,:), 1, [])', c1);
+    x1 = fftfilt(reshape(FIR(2,1,:), 1, [])', c1);
+    y1 = fftfilt(reshape(FIR(3,1,:), 1, [])', c1);
+    z1 = fftfilt(reshape(FIR(4,1,:), 1, [])', c1);
+    c2 = wavread(["../cal/2016.04.20/export/a-format-0" sprintf("%02.0f", m) "-" char(capnames(2)) ".wav"])(range);
+    w2 = fftfilt(reshape(FIR(1,2,:), 1, [])', c2);
+    x2 = fftfilt(reshape(FIR(2,2,:), 1, [])', c2);
+    y2 = fftfilt(reshape(FIR(3,2,:), 1, [])', c2);
+    z2 = fftfilt(reshape(FIR(4,2,:), 1, [])', c2);
+    c3 = wavread(["../cal/2016.04.20/export/a-format-0" sprintf("%02.0f", m) "-" char(capnames(3)) ".wav"])(range);
+    w3 = fftfilt(reshape(FIR(1,3,:), 1, [])', c3);
+    x3 = fftfilt(reshape(FIR(2,3,:), 1, [])', c3);
+    y3 = fftfilt(reshape(FIR(3,3,:), 1, [])', c3);
+    z3 = fftfilt(reshape(FIR(4,3,:), 1, [])', c3);
+    c4 = wavread(["../cal/2016.04.20/export/a-format-0" sprintf("%02.0f", m) "-" char(capnames(4)) ".wav"])(range);
+    w4 = fftfilt(reshape(FIR(1,4,:), 1, [])', c4);
+    x4 = fftfilt(reshape(FIR(2,4,:), 1, [])', c4);
+    y4 = fftfilt(reshape(FIR(3,4,:), 1, [])', c4);
+    z4 = fftfilt(reshape(FIR(4,4,:), 1, [])', c4);
+
+    if (raw == 1)
+      wc = c1 + c2 + c3 + c4;
+      xc = c1 + c2 - c3 - c4;
+      yc = c1 - c2 + c3 - c4;
+      zc = c1 - c2 - c3 + c4;
+      rawlabel = "_raw";
+    else
+      wc = w1 + w2 + w3 + w4;
+      xc = x1 + x2 + x3 + x4;
+      yc = y1 + y2 + y3 + y4;
+      zc = z1 + z2 + z3 + z4;
+      rawlabel = "";
+    end
+    
+    f = figure;
+    hold on;
+    Nfft = 4096;
+    wfft = 20*log10(abs(fft(wc, Nfft)));
+    xfft = 20*log10(abs(fft(xc, Nfft)));
+    yfft = 20*log10(abs(fft(yc, Nfft)));
+    zfft = 20*log10(abs(fft(zc, Nfft)));
+    smooth = 12;
+    wcr = fastsmooth(wfft, smooth, 2, 1);
+    xcr = fastsmooth(xfft, smooth, 2, 1);
+    ycr = fastsmooth(yfft, smooth, 2, 1);
+    zcr = fastsmooth(zfft, smooth, 2, 1);
+    fr = (1:Nfft/2) * Fs/Nfft;
+    grid on;
+    angle = 360 / 16 * (m -1); 
+    title(sprintf("azimuth = %05.1f degrees", angle));
+    axis([300, 20000, -60, 6]);
+    semilogx(fr, wcr(1:end/2), 'k;W;');
+    semilogx(fr, xcr(1:end/2), 'r;X;');
+    semilogx(fr, ycr(1:end/2), 'b;Y;');
+    semilogx(fr, zcr(1:end/2), 'g;Z;');
+    l = legend('W', 'X', 'Y', 'Z');
+    legend(l, 'location', 'southwest');
+    hold off;
+    print(f, '-dpng','-color',["figures/sphear_b_vs_angle_" sprintf("%03.0f", angle) rawlabel ".png"]);
+  end
+
+end
+