remove unused plotting code, scale all plots so that 0 dB is relative to W signal power at a reference frequency range, add optional title argument

calibration/src/sphear_bf_plot_frequency_response.m

@@ -5,18 +5,17 @@
 %% sphear_bf_plot_frequency_response: plot frequency response of B-Format signals
 %%
 %% usage:
-%%   [] = sphear_bf_plot_frequency_response(BF, fmin, fmax, ymin, ymax, Fs, Nfft, plot_on, smooth)
+%%   [] = sphear_bf_plot_frequency_response(BF, frefmin, frefmax, fmin, fmax, ymin, ymax, FS, fft_size, smooth)
 %%
 %% BF: B-Format signals for each measurement
+%% frefmin, frefmax: reference frequency range to normalize W to 0dB
 %% fmin, fmax: range of frequencies to use
 %% ymin, ymax: y range to plot
-%% Fs: sampling rate
-%% Nfft: size of fft used to calculate frequency response
+%% FS: sampling rate
+%% fft_size: size of fft used to calculate frequency 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-2017, Fernando Lopez-Lezcano, All Rights Reserved
@@ -35,82 +34,54 @@
 %%  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_bf_plot_frequency_response(BF, fmin, fmax, ymin, ymax, Fs, Nfft, plot_on, smooth)
-
-  %% no smoothing or plotting by default
+function [] = sphear_bf_plot_frequency_response(BF, frefmin, frefmax, fmin, fmax, ymin, ymax, FS, fft_size, figure_title, smooth)
+  %% no smoothing 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
-
+  %% find average W power in reference frequency range
+  SP = sphear_signal_power(BF, FS, frefmin, frefmax, fft_size);
+  W_ref = mean(SP(:,1));
+  %% calculate B format frequency responses
   for m = 1:size(BF, 1)
-    %% calculate B format frequency responses
     for b = 1:size(BF, 2)
       if (smooth > 0)
-	bfft(m,b,:) = fastsmooth(20*log10(abs(fft(BF(m,b,:), Nfft))), smooth, 2, 1);
+	bfft(m,b,:) = fastsmooth(20*log10(abs(fft(BF(m,b,:) ./ W_ref, fft_size))), smooth, 2, 1);
       else
-	bfft(m,b,:) = 20*log10(abs(fft(BF(m,b,:), Nfft)));
+	bfft(m,b,:) = 20*log10(abs(fft(BF(m,b,:) ./ W_ref, fft_size)));
       end
     end
-
-    %% 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);
-      %% we add the title as a subfigure in the paper
-      %% title(sprintf("azimuth = %.1f degrees", angle));
-      axis([fmin, fmax, ymin, ymax]);
-      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]);
-      ca=get (gcf, "currentaxes");
-      set(ca,"fontweight","bold","linewidth",1.5)
-      pbaspect([16, 9, 1])
-      semilogx(fr, bfft(m,1,1:end/2), wopts, "linewidth", 1.5);
-      semilogx(fr, bfft(m,2,1:end/2), xopts, "linewidth", 1.5);
-      semilogx(fr, bfft(m,3,1:end/2), yopts, "linewidth", 1.5);
-      semilogx(fr, bfft(m,4,1:end/2), zopts, "linewidth", 1.5);
-      if (size(BF, 2) > 4)
-	semilogx(fr, bfft(m,5,1:end/2), 'r--;U;', "linewidth", 1.5);
-	semilogx(fr, bfft(m,6,1:end/2), 'b--;V;', "linewidth", 1.5);
-	l = legend('W', 'X', 'Y', 'Z', 'U', 'V');
-      else
-	l = legend('W', 'X', 'Y', 'Z');
-      end
-      legend(l, 'location', 'west');
-      hold off;
-      if (strcmp(plot_on, 'eps') || strcmp(plot_on, 'epsc') || plot_on == 1)
-	print(f, '-deps', '-tight', copt, ["../data/current/figures/sphear_bf_azimuth_" sprintf("%03.0f", angle) cext ".eps"]);
-      else
-	print(f, '-dpng', copt, ["../data/current/figures/sphear_bf_azimuth_" sprintf("%03.0f", angle) cext ".png"]);
-      end
+  end
+  %% plot frequency responses for each measurement
+  for m = 1:size(BF, 1)
+    fig = figure;
+    hold on;
+    fr = (1:fft_size/2) * FS/fft_size;
+    grid on;
+    angle = 360 / 16 * (m - 1);
+    if (exist("figure_title"))
+      title([figure_title sprintf(" azimuth = %.1f degrees", angle)]);
+    end
+    axis([fmin, fmax, ymin, ymax]);
+    tics('y', [9, 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]);
+    ca=get(gcf, "currentaxes");
+    set(ca, "fontweight", "bold", "linewidth", 1.5);
+    pbaspect([16, 9, 1]);
+    semilogx(fr, bfft(m,1,1:end/2), 'k;W;', "linewidth", 1.5);
+    semilogx(fr, bfft(m,2,1:end/2), 'r;X;', "linewidth", 1.5);
+    semilogx(fr, bfft(m,3,1:end/2), 'b;Y;', "linewidth", 1.5);
+    semilogx(fr, bfft(m,4,1:end/2), 'g;Z;', "linewidth", 1.5);
+    if (size(BF, 2) > 4)
+      semilogx(fr, bfft(m,5,1:end/2), 'r--;U;', "linewidth", 1.5);
+      semilogx(fr, bfft(m,6,1:end/2), 'b--;V;', "linewidth", 1.5);
+      l = legend('W', 'X', 'Y', 'Z', 'U', 'V');
+    else
+      l = legend('W', 'X', 'Y', 'Z');
     end
+    legend(l, 'location', 'west');
+    hold off;
+    %% write file with eps plot
+    print(fig, '-deps', '-tight', '-color', ["../data/current/figures/sphear_bf_azimuth_" sprintf("%03.0f", angle) ".eps"]);
   end
-
 end
-