vectorize sphear_fir_render, add smoothing and size of FFT as parameters, write four channel b format soundfile

sphear_fir_render.m

@@ -3,13 +3,15 @@
 %% sphear_fir_render: render impulse responses throught the FIR filter A2B matrix
 %%
 %% usage:
-%%   [] = sphear_fir_render(FIR, fmin, fmax, N, Npre)
+%%   [] = sphear_fir_render(FIR, fmin, fmax, N, Npre, Nfft, Fs, smooth, plt)
 %%
 %% FIR: A2B matrix of FIR filters
 %% fmin, fmax: range of frequencies to use
 %% N: number of samples to filter
 %% Npre: number of samples before the main impulse peak
+%% Nfft: size of fft used to calculate frequency response
 %% Fs: sampling rate
+%% smooth: number of samples to average
 %% plt: plot response
 %%
 %% Copyright 2016, Fernando Lopez-Lezcano, All Rights Reserved
@@ -28,7 +30,7 @@
 %%  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 [BFF] = sphear_fir_render(FIR, fmin, fmax, N, Npre, Fs, plt)
+function [BFF] = sphear_fir_render(FIR, fmin, fmax, N, Npre, Nfft, Fs, smooth, plt)
 
   [ref, Fs] = wavread('../cal/2016.04.20/export/reference-001-1.wav');
   %% find highest peak location and define sample range
@@ -41,121 +43,67 @@ function [BFF] = sphear_fir_render(FIR, fmin, fmax, N, Npre, Fs, plt)
   
   capnames = cellstr(['LF'; 'RF'; 'LB'; 'RB']);
 
-  raw = 0;
   writefile = 1;
 
-  minindex = N;
-  maxpeak = 0;
+  idx = 1;
   
   for m = 1:16
-    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);
+    %% render b format components
+    bf(idx,:,:) = sphear_a_to_b(m, FIR, range, N);
 
-    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);
-
-    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);
-
-    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);
-      rawlabel = "_raw";
-    else
-      wc = (w1 + w2 + w3 + w4);
-      xc = (x1 + x2 + x3 + x4);
-      yc = (y1 + y2 + y3 + y4);
-      zc = (z1 + z2 + z3 + z4);
-      rawlabel = "";
+    %% calculate B format frequency responses
+    for b = 1:4
+      bfft(idx,b,:) = fastsmooth(20*log10(abs(fft(bf(idx,b,:), Nfft))), smooth, 2, 1);
     end
 
-    BFF(:,1,m) = wc;
-    BFF(:,2,m) = xc;
-    BFF(:,3,m) = yc;
-    BFF(:,4,m) = zc;
-    
+    %% write soundfile with rendered B format signals
     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"]);
+      %% if we write a 32 bit soundfile we need to scale by
+      %% 256 if want the same amplitude as when writing a
+      %% 16 bit sample soundfile...
+      wavwrite(reshape(bf(idx,:,:),4,[])' / 256, Fs, 32, ["../cal/2016.04.20/export/b-format-" 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([fmin, fmax, -45, 0]);
+      axis([fmin, fmax, -45, 12]);
       tics('y', [6, 0, -2.5, -5, -7.5, -10, -12.5, -15, -17.5, -20, -22.5, -25, -30, -40, -50]);
       tics('x', [300, 500, 1000, 2000, 3000, 5000, 7000, 10000, 15000, 20000]);
-      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;');
+      semilogx(fr, bfft(idx,1,1:end/2), 'k;W;');
+      semilogx(fr, bfft(idx,2,1:end/2), 'r;X;');
+      semilogx(fr, bfft(idx,3,1:end/2), 'b;Y;');
+      semilogx(fr, bfft(idx,4,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"]);
+      print(f, '-deps','-color',["figures/sphear_b_vs_angle_" sprintf("%03.0f", angle) ".eps"]);
     end
+
+    idx = idx + 1;
   end
 
-  [minindex, maxpeak]
+end
+
+
+function [bf] = sphear_a_to_b(m, FIR, range, N);
+  capnames = cellstr(['LF'; 'RF'; 'LB'; 'RB']);
+  %% read in capsule signals and apply FIR filter
+  for c = 1:4
+    cs(c,:) = wavread(["../cal/2016.04.20/export/a-format-0" sprintf("%02.0f", m) "-" char(capnames(c)) ".wav"]);
+    cr(c,:) = tukeywin(size(range)(2), 0.0) .* cs(c,:)'(range);
+    for b = 1:4
+      bs(c,b,:) = fftfilt(reshape(FIR(b,c,:), 1, [])', cr(c,:), N);
+    end
+  end
+  %% add filtered capsule output to b format components
+  for b = 1:4
+    bf(b,:) = bs(1,b,:) + bs(2,b,:) + bs(3,b,:) + bs(4,b,:);
+  end
 end