%%
%% *SpHEAR Project
%% sphear_fir_render: render impulse responses throught the FIR filter A2B matrix
%%
%% usage:
%%   [] = sphear_fir_render(FIR, fmin, fmax, N, Npre, Nfft, Fs, smooth, plt)
%%
%% 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
%% Nfft: size of fft used to calculate frequency response
%% Fs: sampling rate
%% smooth: number of samples to average
%% plt: plot response
%%
%% 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(FIR, fmin, fmax, N, Npre, Nfft, Fs, smooth, plt)

  [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(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']);

  writefile = 1;

  idx = 1;
  
  for m = 1:16
    %% 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

    %% write soundfile with rendered B format signals
    if (writefile == 1)
      %% if we write a 32 bit soundfile we need to scale by
      %% 256 if want the same amplitude as when writing a
      %% 16 bit sample soundfile...
      wavwrite(reshape(bf(idx,:,:),4,[])' / 256, Fs, 32, ["../cal/2016.04.20/export/b-format-" sprintf("%02.0f", m) ".wav"]);
    end

    if (plt == 1)
      f = figure;
      hold on;
      fr = (1:Nfft/2) * Fs/Nfft;
      grid on;
      angle = 360 / 16 * (m - 1); 
      title(sprintf("azimuth = %05.1f degrees", angle));
      axis([fmin, fmax, -45, 12]);
      tics('y', [6, 0, -2.5, -5, -7.5, -10, -12.5, -15, -17.5, -20, -22.5, -25, -30, -40, -50]);
      tics('x', [300, 500, 1000, 2000, 3000, 5000, 7000, 10000, 15000, 20000]);
      semilogx(fr, bfft(idx,1,1:end/2), 'k;W;');
      semilogx(fr, bfft(idx,2,1:end/2), 'r;X;');
      semilogx(fr, bfft(idx,3,1:end/2), 'b;Y;');
      semilogx(fr, bfft(idx,4,1:end/2), 'g;Z;');
      l = legend('W', 'X', 'Y', 'Z');
      legend(l, 'location', 'southwest');
      hold off;
      print(f, '-deps','-color',["figures/sphear_b_vs_angle_" sprintf("%03.0f", angle) ".eps"]);
    end

    idx = idx + 1;
  end

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