add code to plot the capsule polar response as a function of frequency

capsule_plot.m

@@ -0,0 +1,58 @@
+%%
+%% *SpHEAR Project
+%%
+%% capsule_plot: plot the polar response of a capsule
+%%
+%% usage:
+%%   capsule_plot(R, fmin, fmax)
+%%
+%% R: capsule response in frequency range
+%% 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 capsule_plot(i, R, fmin, fmax, rmax, col)
+
+  m = size(R, 1);
+  n = 32;
+  theta = (0:m-1)*(pi/(m-1));
+  thetaf = (0:n)*(pi/n);
+  
+  hold on;
+  axis([-rmax, rmax, 0, rmax], "equal", "tickx");
+  
+  title(["EM182 polar pattern"]);
+
+  if (i == 1)
+    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
+  end
+  
+  [x, y] = pol2cart(theta', abs(R(:)));
+  p = plot(x, y, [".-;" sprintf("%d to %d Hz", fmin, fmax) ";"], 'LineWidth', 0.5);
+  set(p, 'color', col); 
+
+end
+

capsule_polar_response.m

@@ -0,0 +1,62 @@
+%%
+%% *SpHEAR Project
+%%
+%% capsule_polar_response: plot polar response of capsule for different frequency bands
+%%
+%% usage:
+%%   [] = capsule_polar_respones()
+%%
+%% The names of reference microphone impulse response and capsule measurements are
+%% currently hardwired into this file, this needs to be fixed and turned into
+%% input parameters
+%%
+%% 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 [CAP] = capsule_polar_response(measurements, fmin, fmax, N, Npre)
+  pkg load signal;
+
+  %% read in reference microphone, select samples to process
+  [ref, Fs] = wavread('../cal/2016.04.20/export/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), N);
+  
+  %% indexes in array for frequency bins
+  imin = round((N/2) * fmin / (Fs/2));
+  imax = round((N/2) * fmax / (Fs/2));
+
+  capnames = cellstr(['LF'; 'RF'; 'LB'; 'RB']);
+
+  %% all submitted measurements
+  i = 1;
+  for m = measurements
+    %% read in capsule measurement and equalize to reference microphone
+    %% this is only valid for the frequency range covered by the test speaker,
+    %% anything outside is undefined
+    %%
+    capsule = wavread(["../cal/2016.04.20/export/capsule-0" sprintf("%02.0f", m) "-2.wav"]);
+    equalized = fft(capsule(range), N) ./ ref_fft;
+    
+    %% return total power of signal in frequency range
+    CAP(i,:) = sqrt(sum(equalized(imin:imax).*conj(equalized(imin:imax)))/size(equalized(imin:imax),1));
+    i = i + 1;
+  end
+
+end
+

capsule_polar_response_all.m

@@ -0,0 +1,61 @@
+%%
+%% *SpHEAR Project
+%%
+%% capsule_polar_response_all: calculate polar response of capsule for frequency ranges
+%%
+%% usage:
+%%   [] = capsule_polar_response_all(m, fmin, fmax, Noct, Fs, Nsamp, Npre, plt)
+%%
+%% m: vector of measurement numbers to take into account
+%% 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
+%% Fs: sampling rate
+%% Nsamp: time domain samples to use
+%% Npre: samples to use before the main impulse peak
+%%
+%% 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 [] = capsule_polar_response_all(m, fmin, fmax, Noct, Fs, Nsamp, Npre, plt)
+
+  %% 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)';
+
+  colors = 'brgkmc';
+  f = figure;
+  title(["EM182 capsule polar pattern"]);
+
+  for i = (1:size(fc))
+    R = capsule_polar_response(m, fe(i), fe(i + 1), Nsamp, Npre, 1.5);
+    if plt
+      capsule_plot(1, R, fe(i), fe(i + 1), 1.8, colors(mod(i - 1, size(colors)(2)) + 1));
+    end
+  end
+  legend('location', 'west');
+  print(f, '-depsc', 'figures/capsule_polar_response.eps');
+end
+