add functions to plot polar response of a B format signal

sphear_wxyz_plot.m

@@ -0,0 +1,92 @@
+%%
+%% *SpHEAR Project
+%%
+%% sphear_wxyz_plot: plot the polar response of a B format signal in a frequency range
+%%
+%% usage:
+%%   sphear_bf_plot(BF, fmin, fmax)
+%%
+%% BF: B format signal array
+%% fmin, fmax: range of frequencies to consider (to title plot)
+%%
+%% Copyright Fernando Lopez-Lezcano 2016
+%% 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_wxyz_plot(BF, fmin, fmax, rmax)
+
+  if (!exist("plot_on"))
+    plot_on = "epsc";
+  else
+    if (plot_on == 1)
+      plot_on = "epsc"
+    end
+  end
+
+  %% to close the polar plots we repeat the first point
+  m = size(BF, 1);
+  n = 32;
+  BF(m + 1,:) = BF(1,:);
+  theta = (0:m)*(2*pi/m);
+  thetaf = (0:n)*(2*pi/n);
+  
+  %% plot it as a function of angle
+  bname = cellstr(['W'; 'X'; 'Y'; 'Z']);
+
+  f = figure;
+  for b = 1:4
+    subplot(2, 2, b);
+    hold on;
+    grid off;
+    axis square;
+    axis([-rmax, rmax, -rmax, rmax]);
+    ca=get (gcf, "currentaxes");
+    set(ca,"fontweight","bold","linewidth",1.5)
+    title([char(bname(b)) ": " sprintf("%05.0d", fmin) " to " sprintf("%05.0d", fmax) " Hz"]);
+
+    if (strcmp(plot_on, "epsc") || strcmp(plot_on, "pngc"))
+      %% color plot
+      dopts = '.-r';
+      ropts = '-b';
+      copt = '-color';
+      cext = "_color";
+    else
+      %% mono plot
+      dopts = '.-';
+      ropts = '-';
+      copt = '-mono';
+      cext = "_mono";
+    end
+
+    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);
+    end
+    %% radial lines
+    for ang = theta
+      plot([0 rmax*cos(ang)], [0 rmax*sin(ang)], ':', 'Color', grey);
+    end
+
+    [xr, yr] = pol2cart(thetaf', ones(n + 1,1));
+    plot(xr, yr, ropts);
+    [x, y] = pol2cart(theta', abs(BF(:,b)));
+    plot(x, y, dopts, 'linewidth', 1.5, 'markersize', 6);
+    
+  end
+  print(f, '-deps','-tight',copt,["figures/sphear_polar_multi_" sprintf("%05d", fmin) "_" sprintf("%05d", fmax) cext ".eps"]);
+end
+

sphear_wxyz_plot_all.m

@@ -0,0 +1,63 @@
+%%
+%% *SpHEAR Project
+%%
+%% sphear_wxyz_plot_all: plot the polar response of a B format signal
+%%
+%% usage:
+%%   [B] = sphear_wxyz_plot_all(fmin, fmax, Noct, Nbuf, Npre, Nfft, plt)
+%%
+%% 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
+%% Nbuf: total number of samples to analyze
+%% Npre: number of samples before the main impulse
+%% Nfft: size of fft
+%% plt: 0 -> just return the measurements, 1 -> do plots
+%%
+%% 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 [B] = sphear_wxyz_plot_all(fmin, fmax, Noct, Nbuf, Npre, Nfft, plt)
+
+  plotz = 0;
+  measurements = (1:16);
+  
+  %% 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_wxyz_power(measurements, 600, 1200, Nbuf, Npre, Nfft, 1);
+  norm = mean(ref(:,1));
+
+  for i = (1:size(fc))
+    B(:,:,i) = sphear_wxyz_power(measurements, fe(i), fe(i + 1), Nbuf, Npre, Nfft, plotz) / norm;
+    if (plt == 1)
+      sphear_wxyz_plot(B(:,:,i), fe(i), fe(i + 1), 1.5);
+    end
+  end
+
+end
+

sphear_wxyz_power.m

@@ -0,0 +1,79 @@
+%%
+%% *SpHEAR Project
+%%
+%% sphear_a_power: return array of power for WXYZ B format signals in a range of frequencies
+%%
+%% usage:
+%%   [B Fs] = sphear_wxyz_power(measurements, fmin, fmax, N, Npre, Nfft, plotz)
+%%
+%% measurements: measurements to consider
+%% fmin, fmax: range of frequencies to consider
+%% N: total number of samples to analyze
+%% Npre: number of samples before the main impulse
+%% Nfft: size of fft
+%% plotz: 0 -> plot cardioid, 1 -> plot Z
+%%
+%% The names of reference microphone impulse response and capsule measurements are
+%% currently hardwired into this file, this needs to be fixed and turned into
+%% parameters
+%%
+%% 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 [B Fs] = sphear_wxyz_power(measurements, fmin, fmax, N, Npre, Nfft, plotz)
+
+  %% read in reference microphone, select samples to process
+  [ref, Fs] = wavread('../data/2016.08.03/st250/reference-001-1.wav');
+  [peaks, locs] = findpeaks(abs(ref));
+  [~, index] = max(peaks);
+  start = locs(index) - Npre;
+  range = start:(start + N - 1);
+  ref_fft = fft(ref(range), Nfft);
+  
+  %% indexes in array for frequency bins
+  imin = round((N/2) * fmin / (Fs/2));
+  imax = round((N/2) * fmax / (Fs/2));
+
+  %% name of B format components
+  bnames = cellstr(['W'; 'X'; 'Y'; 'Z']);
+
+  %% all measurements
+  for m = measurements
+    %% all components
+    for c = 1:4
+      %% all b format components, 4 channel soundfile
+      bsignal = wavread(["../data/2016.08.03/st250/b-format-0" sprintf("%02.0f", m)  "-" char(bnames(c)) ".wav"]);
+      if (c == 1)
+	%% FUMA to SN3D for W
+	bsignal = bsignal / 0.707;
+      end
+      bfft = fft(bsignal(range), Nfft) ./ ref_fft;
+      if (plotz == 0)
+	%% calculate a forward facing cardiod
+	cfft = (fft(bsignal(range) + 1.0 * bsignal(range), Nfft)) ./ ref_fft;
+      end
+      %% return total power of signal in frequency range
+      B(m, c) = sqrt(sum(bfft(imin:imax,:).*conj(bfft(imin:imax,:)))/size(bfft(imin:imax,:),1));
+      %% make a cardiod if we don't want Z
+      if (c == 4)
+	if (plotz == 0)
+	  B(m, 4) = sqrt(sum(cfft(imin:imax).*conj(cfft(imin:imax)))/size(cfft(imin:imax),1)) / 2.0;
+	end
+      end
+    end
+  end
+end
+