add polar plot responses with arbitrary az/el angles

calibration/src/sphear_bf_az_el_plot.m

@@ -0,0 +1,128 @@
+%% -*- mode: octave -*-
+%%
+%% *SpHEAR Project
+%%
+%% sphear_bf_plot: plot the polar response of a B format signal
+%%
+%% usage:
+%%   [BF] = sphear_bf_plot(BF, AZ, EL, fmin, fmax)
+%%
+%% BF: B format signal array
+%% M: projection matrix in measurement directions
+%% AZ: azimuth angles
+%% EL: elevation angles
+%% fmin, fmax: range of frequencies to consider (to title plot)
+%% cardioid: whether we have Z or a cardioid in its place
+%% subplots: whether we want individual plots or subplots
+%%
+%% Copyright Fernando Lopez-Lezcano 2016-2018
+%% 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 [BF] = sphear_bf_az_el_plot(BF, M, AZ, EL, fmin, fmax, rmax, cardioid = 0, subplots = 1)
+  %% sort azimuth angles
+  [AZS, I] = sort(AZ);
+  ELS = EL(I);
+  BFS(:,:) = BF(I,:);
+  MS(:,:) = M(I,:);
+  %% to close the polar plots we repeat the first point
+  m = size(BF, 1);
+  BFS(m + 1,:) = BFS(1,:);
+  MS(m + 1,:) = MS(1,:);
+  AZS(m + 1) = AZS(1);
+  ELS(m + 1) = ELS(1);
+  %% angles for plotting polar grid
+  n = 32;
+  theta = (0:m)*(2*pi/m);
+  thetaf = (0:n)*(2*pi/n);
+  %% plot it as a function of angle
+  if (cardioid)
+    bname = cellstr(['W'; 'X'; 'Y'; 'Cardioid'; 'U'; 'V']);
+  else
+    if (size(BF, 2) > 6)
+      bname = cellstr(['W'; 'X'; 'Y'; 'Z'; 'S'; 'T'; 'U'; 'V']);
+    else
+      bname = cellstr(['W'; 'X'; 'Y'; 'Z'; 'U'; 'V']);
+    end
+  end
+  %% go and plot
+  f = figure;
+  for b = 1:size(BF, 2)
+    if (subplots == 0)
+      f = figure;
+      font_size = 28;
+      line_w = 1.0;
+      line_W = 3;
+      dot_size = 15;
+      title([char(bname(b))]);
+    else
+      if (size(BF,2) > 4)
+	subplot(2, 4, b);
+      else
+	subplot(2, 2, b);
+      end
+      font_size = 18;
+      line_w = 0.5;
+      line_W = 1.5;
+      dot_size = 10;
+      title([char(bname(b)) sprintf(" (%.2f - %.2f Hz)", fmin, fmax)]);
+    end
+    hold on;
+    grid off;
+    axis("off", "equal", "nolabel");
+    ca=get (gcf, "currentaxes");
+    set(ca,"fontweight", "bold", "fontsize", font_size)
+    grey = [0.1, 0.1, 0.1];
+    %% radius circles
+    for rad = 0.25:0.25:rmax
+      [xr, yr] = pol2cart(thetaf', ones(n + 1,1) * rad);
+      plot(xr, yr, '--', 'color', grey, 'linewidth', line_w);
+    end
+    %% radial lines
+    for ang = theta
+      plot([0 rmax*cos(ang)], [0 rmax*sin(ang)], '--', 'color', grey, 'linewidth', line_w);
+    end
+    [xr, yr] = pol2cart(thetaf', ones(n + 1,1));
+    plot(xr, yr, "b-", "linewidth", line_w);
+    %% plot real harmonic normalized to max measured harmonic
+    if (b > 0)
+      if (b == 1)
+	%% M does not have W
+	db_e = 1.0;
+      else
+	%% real harmonic scaled by maximum of measured harmonic
+	for m = 1:size(MS,1)
+	  if (max(abs(MS(:,b-1))) > 0)
+	    scaled = 0;
+	    if (scaled == 1)
+	      db_e(m,:) = abs(MS(m,b-1)) * (max(abs(BFS(:,b))) / max(abs(MS(:,b-1))));
+	    else
+	      db_e(m,:) = abs(MS(m,b-1));
+	    end
+	  else
+	    db_e(m,:) = 0.0;
+	  end
+	end
+      end
+      db_e_max = max(abs(db_e));
+      [xm, ym] = pol2cart(AZS, abs(db_e));
+      plot(-ym, xm, "k.--", 'linewidth', line_W, 'markersize', dot_size);
+    end
+    %% plot measured harmonic
+    [x, y] = pol2cart(AZS, abs(BFS(:,b)));
+    plot(-y, x, "r.-", 'linewidth', line_W, 'markersize', dot_size);
+  end
+end
+

calibration/src/sphear_bf_az_el_plot_polar_response.m

@@ -0,0 +1,79 @@
+%% -*- mode: octave -*-
+%%
+%% *SpHEAR Project
+%%
+%% sphear_bf_az_el_plot_polar_response: plot polar response of B format components
+%%
+%% usage:
+%%   sphear_bf_plot_polar_response(B, M, AZ, EL, frefmin, frefmax, fmin, fmax, Noct, Nfft, Fs, band, cardioid = 0, azim = 45)
+%%
+%% B: B format signal
+%% M: projections matrix
+%% AZ: azimuth angles
+%% EL: elevation angles
+%% frefmin, frefmax: reference frequencies for power reference
+%% fmin, fmax: lower and upper frequency limits
+%% Noct: width of frequency ranges in fractions of an octave, this
+%%       will be changed slightly so that the frequency range is
+%%       covered by equally spaced frequency bands
+%% Nfft: FFT size
+%% Fs: sampling rate
+%% band: single or array of frequencies, they will be matched to the bands
+%%
+%% optional:
+%%
+%% cardioid: replace Z with a synthetic cardioid
+%% azimuth: azimuth angle for the synthetic cardioid
+%%
+%% Copyright 2016-2018, 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_bf_az_el_plot_polar_response(B, M, AZ, EL, frefmin, frefmax, fmin, fmax, Noct, Nfft, Fs, band, cardioid = 0, azim = 45, subplots = 1)
+  %% calculate number of frequency bands spaced Noct
+  n = floor(log(fmax/fmin)/log(2^(1/Noct)));
+  %% adjust width of band so they fit exactly in frequency range
+  N = log(2)/(log(fmax/fmin)/n);
+  printf("%d frequency bands, %f bands per octave\n", n, N);
+
+  %% first center frequency
+  fcmin = fmin * 2^(1/(N * 2));
+  %% all center and edge frequencies
+  fc = fcmin * 2.^((0:n - 1)/N)';
+  fe = fmin * 2.^((0:n)/N)';
+
+  %% normalize amplitude with respect to W between 600 and 1200 Hz
+  ref = sphear_signal_power(B, Fs, frefmin, frefmax, Nfft);
+  norm = mean(ref(:,1));
+  
+  %% plot all bands if we receive an empty array
+  if (numel(band) == 0)
+    band = fc;
+  end
+
+  idx = 1;
+  for ib = 1:numel(band)
+    printf("processing band around %d\n", band(ib));
+    for i = (1:size(fc))
+      if ((fe(i) <= band(ib)) && ((fe(i+1) >= band(ib))))
+	printf("fc = %d (between %d and %d)\n", fc(i), fe(i), fe(i+1));
+	BP = sphear_signal_power(B, Fs, fe(i), fe(i + 1), Nfft, cardioid, azim) / norm;
+	sphear_bf_az_el_plot(BP, M, AZ, EL, fe(i), fe(i + 1), 1.25, cardioid, subplots);
+      end
+      idx = idx + 1;
+    end
+  end
+end
+