sphear_fir_render: many changes, write b format files, find and report max signal for components

sphear_fir_render.m

@@ -9,6 +9,8 @@
 %% fmin, fmax: range of frequencies to use
 %% N: number of samples to filter
 %% Npre: number of samples before the main impulse peak
+%% Fs: sampling rate
+%% plt: plot response
 %%
 %% Copyright 2016, Fernando Lopez-Lezcano, All Rights Reserved
 %% nando@ccrma.stanford.edu
@@ -26,82 +28,139 @@
 %%  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_fir_render(FIR, fmin, fmax, N, Npre)
+function [BFF] = sphear_fir_render(FIR, fmin, fmax, N, Npre, Fs, plt)
 
   [ref, Fs] = wavread('../cal/2016.04.20/export/reference-001-1.wav');
   %% find highest peak location and define sample range
-  [peaks, locs] = findpeaks(abs(ref));
-  [~, index] = max(peaks);
-  start = locs(index) - Npre;
-  range = start:(start + N - 1);
+  [index, value] = findmaxpeak(ref);
+  r_start = index - Npre;
+  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']);
 
   raw = 0;
+  writefile = 1;
 
+  minindex = N;
+  maxpeak = 0;
+  
   for m = 1:16
-    c1 = wavread(["../cal/2016.04.20/export/a-format-0" sprintf("%02.0f", m) "-" char(capnames(1)) ".wav"])(range);
+    c1r = wavread(["../cal/2016.04.20/export/a-format-0" sprintf("%02.0f", m) "-" char(capnames(1)) ".wav"]);
+    c1 = c1r(range);
+    [index, peak] = findmaxpeak(c1r);
+    maxpeak = max(maxpeak, abs(peak));
+    minindex = min(minindex, (index - r_start));
     w1 = fftfilt(reshape(FIR(1,1,:), 1, [])', c1);
     x1 = fftfilt(reshape(FIR(2,1,:), 1, [])', c1);
     y1 = fftfilt(reshape(FIR(3,1,:), 1, [])', c1);
     z1 = fftfilt(reshape(FIR(4,1,:), 1, [])', c1);
-    c2 = wavread(["../cal/2016.04.20/export/a-format-0" sprintf("%02.0f", m) "-" char(capnames(2)) ".wav"])(range);
+
+    c2r = wavread(["../cal/2016.04.20/export/a-format-0" sprintf("%02.0f", m) "-" char(capnames(2)) ".wav"]);
+    c2 = c2r(range);
+    [index, peak] = findmaxpeak(c2r);
+    maxpeak = max(maxpeak, peak);
+    minindex = min(minindex, (index - r_start));
     w2 = fftfilt(reshape(FIR(1,2,:), 1, [])', c2);
     x2 = fftfilt(reshape(FIR(2,2,:), 1, [])', c2);
     y2 = fftfilt(reshape(FIR(3,2,:), 1, [])', c2);
     z2 = fftfilt(reshape(FIR(4,2,:), 1, [])', c2);
-    c3 = wavread(["../cal/2016.04.20/export/a-format-0" sprintf("%02.0f", m) "-" char(capnames(3)) ".wav"])(range);
+
+    c3r = wavread(["../cal/2016.04.20/export/a-format-0" sprintf("%02.0f", m) "-" char(capnames(3)) ".wav"]);
+    c3 = c3r(range);
+    [index, peak] = findmaxpeak(c3r);
+    maxpeak = max(maxpeak, peak);
+    minindex = min(minindex, (index - r_start));
     w3 = fftfilt(reshape(FIR(1,3,:), 1, [])', c3);
     x3 = fftfilt(reshape(FIR(2,3,:), 1, [])', c3);
     y3 = fftfilt(reshape(FIR(3,3,:), 1, [])', c3);
     z3 = fftfilt(reshape(FIR(4,3,:), 1, [])', c3);
-    c4 = wavread(["../cal/2016.04.20/export/a-format-0" sprintf("%02.0f", m) "-" char(capnames(4)) ".wav"])(range);
+
+    c4r = wavread(["../cal/2016.04.20/export/a-format-0" sprintf("%02.0f", m) "-" char(capnames(4)) ".wav"]);
+    c4 = c4r(range);
+    [index, peak] = findmaxpeak(c4r);
+    maxpeak = max(maxpeak, peak);
+    minindex = min(minindex, (index - r_start));
     w4 = fftfilt(reshape(FIR(1,4,:), 1, [])', c4);
     x4 = fftfilt(reshape(FIR(2,4,:), 1, [])', c4);
     y4 = fftfilt(reshape(FIR(3,4,:), 1, [])', c4);
     z4 = fftfilt(reshape(FIR(4,4,:), 1, [])', c4);
 
     if (raw == 1)
-      wc = c1 + c2 + c3 + c4;
-      xc = c1 + c2 - c3 - c4;
-      yc = c1 - c2 + c3 - c4;
-      zc = c1 - c2 - c3 + c4;
+      wc = (c1 + c2 + c3 + c4);
+      xc = (c1 + c2 - c3 - c4);
+      yc = (c1 - c2 + c3 - c4);
+      zc = (c1 - c2 - c3 + c4);
       rawlabel = "_raw";
     else
-      wc = w1 + w2 + w3 + w4;
-      xc = x1 + x2 + x3 + x4;
-      yc = y1 + y2 + y3 + y4;
-      zc = z1 + z2 + z3 + z4;
+      wc = (w1 + w2 + w3 + w4);
+      xc = (x1 + x2 + x3 + x4);
+      yc = (y1 + y2 + y3 + y4);
+      zc = (z1 + z2 + z3 + z4);
       rawlabel = "";
     end
+
+    BFF(:,1,m) = wc;
+    BFF(:,2,m) = xc;
+    BFF(:,3,m) = yc;
+    BFF(:,4,m) = zc;
     
-    f = figure;
-    hold on;
-    Nfft = 4096;
-    wfft = 20*log10(abs(fft(wc, Nfft)));
-    xfft = 20*log10(abs(fft(xc, Nfft)));
-    yfft = 20*log10(abs(fft(yc, Nfft)));
-    zfft = 20*log10(abs(fft(zc, Nfft)));
-    smooth = 12;
-    wcr = fastsmooth(wfft, smooth, 2, 1);
-    xcr = fastsmooth(xfft, smooth, 2, 1);
-    ycr = fastsmooth(yfft, smooth, 2, 1);
-    zcr = fastsmooth(zfft, smooth, 2, 1);
-    fr = (1:Nfft/2) * Fs/Nfft;
-    grid on;
-    angle = 360 / 16 * (m -1); 
-    title(sprintf("azimuth = %05.1f degrees", angle));
-    axis([300, 20000, -60, 6]);
-    semilogx(fr, wcr(1:end/2), 'k;W;');
-    semilogx(fr, xcr(1:end/2), 'r;X;');
-    semilogx(fr, ycr(1:end/2), 'b;Y;');
-    semilogx(fr, zcr(1:end/2), 'g;Z;');
-    l = legend('W', 'X', 'Y', 'Z');
-    legend(l, 'location', 'southwest');
-    hold off;
-    print(f, '-dpng','-color',["figures/sphear_b_vs_angle_" sprintf("%03.0f", angle) rawlabel ".png"]);
+    if (writefile == 1)
+      maxw = max(abs(wc));
+      maxx = max(abs(xc));
+      maxy = max(abs(yc));
+      maxz = max(abs(zc));
+      maxbf = max([maxw, maxx, maxy, maxz]);
+      printf("m%d: max value is %f\n", m, maxbf);
+      %% don't know why the written soundfile is clipped...
+      %% the max value we calculate is less that 1.0...
+      %% so set up an attenuation factor...
+      att = 0.05;
+      wavwrite(wc * att, Fs, 32, ["../cal/2016.04.20/export/b-format-w-" sprintf("%02.0f", m) ".wav"]);
+      wavwrite(xc * att, Fs, 32, ["../cal/2016.04.20/export/b-format-x-" sprintf("%02.0f", m) ".wav"]);
+      wavwrite(yc * att, Fs, 32, ["../cal/2016.04.20/export/b-format-y-" sprintf("%02.0f", m) ".wav"]);
+      wavwrite(zc * att, Fs, 32, ["../cal/2016.04.20/export/b-format-z-" sprintf("%02.0f", m) ".wav"]);
+    end
+
+    if (plt == 1)
+      f = figure;
+      hold on;
+      Nfft = 4096;
+      wfft = 20*log10(abs(fft(wc, Nfft)));
+      xfft = 20*log10(abs(fft(xc, Nfft)));
+      yfft = 20*log10(abs(fft(yc, Nfft)));
+      zfft = 20*log10(abs(fft(zc, Nfft)));
+      smooth = 12;
+      wcr = fastsmooth(wfft, smooth, 2, 1);
+      xcr = fastsmooth(xfft, smooth, 2, 1);
+      ycr = fastsmooth(yfft, smooth, 2, 1);
+      zcr = fastsmooth(zfft, smooth, 2, 1);
+      fr = (1:Nfft/2) * Fs/Nfft;
+      grid on;
+      angle = 360 / 16 * (m - 1); 
+      title(sprintf("azimuth = %05.1f degrees", angle));
+      axis([300, 20000, -60, 6]);
+      tics('y', [6, 0, -2.5, -5, -7.5, -10, -12.5, -15, -17.5, -20, -22.5, -25, -30, -40, -50]);
+      semilogx(fr, wcr(1:end/2), 'k;W;');
+      semilogx(fr, xcr(1:end/2), 'r;X;');
+      semilogx(fr, ycr(1:end/2), 'b;Y;');
+      semilogx(fr, zcr(1:end/2), 'g;Z;');
+      l = legend('W', 'X', 'Y', 'Z');
+      legend(l, 'location', 'southwest');
+      hold off;
+      print(f, '-dpng','-color',["figures/sphear_b_vs_angle_" sprintf("%03.0f", angle) rawlabel ".png"]);
+    end
   end
 
+  [minindex, maxpeak]
 end
 
+
+function [i, v] = findmaxpeak(signal);
+  [peaks, locs] = findpeaks(signal, "DoubleSided");
+  [~, index] = max(peaks);
+  i = locs(index);
+  v = peaks(index);
+end