optimize parameters, add plot options for eps, png, color and mono, remove functions moved to sphear_common

sphear_fir_render.m

@@ -1,18 +1,22 @@
 %%
 %% *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)
+%%   [] = sphear_fir_render(FIR, fmin, fmax, N, Npre, Nfft, plot_on, 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
 %% Nfft: size of fft used to calculate frequency response
-%% Fs: sampling rate
-%% smooth: number of samples to average
-%% plt: plot response
+%%
+%% 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
 %%
 %% Copyright 2016, Fernando Lopez-Lezcano, All Rights Reserved
 %% nando@ccrma.stanford.edu
@@ -30,7 +34,19 @@
 %%  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)
+function [BFF] = sphear_fir_render(FIR, fmin, fmax, N, Npre, 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
 
   [ref, Fs] = wavread('../cal/2016.04.20/export/reference-001-1.wav');
   %% find highest peak location and define sample range
@@ -39,77 +55,73 @@ function [BFF] = sphear_fir_render(FIR, fmin, fmax, N, Npre, Nfft, Fs, smooth, p
   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;
+  printf("reference: peak %f at %d, sample range [%d %d], sampling rate: %d\n",
+	 value, index, r_start, r_end, Fs);
   
   for m = 1:16
     %% render b format components
-    bf(idx,:,:) = sphear_a_to_b(m, FIR, range, N);
+    bf(m,:,:) = 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);
+      if (smooth > 0)
+	bfft(m,b,:) = fastsmooth(20*log10(abs(fft(bf(m,b,:), Nfft))), smooth, 2, 1);
+      else
+	bfft(m,b,:) = 20*log10(abs(fft(bf(m,b,:), Nfft)));
+      end
     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 we write a 32 bit soundfile we need to scale by 256
+    %% if we want the same amplitude as when writing a 16 bit
+    %% per sample soundfile, otherwise we get nasty clipping
+    wavwrite(reshape(bf(m,:,:),4,[])' / 256, Fs, 32, ["../cal/2016.04.20/export/b-format-" sprintf("%02.0f", m) ".wav"]);
 
-    if (plt == 1)
+    %% 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;';
+	zopts = 'g;Z;';
+	copt = '-color';
+	cext = "_color";
+      else
+	%% mono plot
+	wopts = '-;W;';
+	xopts = '--;X;';
+	yopts = '-.;Y;';
+	zopts = ':;Z;';
+	copt = '-mono';
+	cext = "_mono";
+      end
       f = figure;
       hold on;
       fr = (1:Nfft/2) * Fs/Nfft;
       grid on;
-      angle = 360 / 16 * (m - 1); 
-      title(sprintf("azimuth = %05.1f degrees", angle));
+      angle = 360 / 16 * (m - 1);
+      %% we add the title as a subfigure in the paper
+      %% title(sprintf("azimuth = %.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('y', [6, 3, 0, -3, -6, -9, -12, -15, -20, -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;');
+      pbaspect([16, 9, 1])
+      semilogx(fr, bfft(m,1,1:end/2), wopts);
+      semilogx(fr, bfft(m,2,1:end/2), xopts);
+      semilogx(fr, bfft(m,3,1:end/2), yopts);
+      semilogx(fr, bfft(m,4,1:end/2), zopts);
       l = legend('W', 'X', 'Y', 'Z');
-      legend(l, 'location', 'southwest');
+      legend(l, 'location', 'west');
       hold off;
-      print(f, '-deps','-color',["figures/sphear_b_vs_angle_" sprintf("%03.0f", angle) ".eps"]);
+      if (strcmp(plot_on, 'eps') || strcmp(plot_on, 'epsc') || plot_on == 1)
+	print(f, '-deps', '-tight', copt, ["figures/sphear_bf_azimuth_" sprintf("%03.0f", angle) cext ".eps"]);
+      else
+	print(f, '-dpng', copt, ["figures/sphear_bf_azimuth_" sprintf("%03.0f", angle) cext ".png"]);
+      end
     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