sphear_fir_design: add tentative code to add an optional high frequency boost; create running averages for the A2B matrices and WXYZ filter shapes - both going up to the transition frequency and also down so that the transition between the two modes is smoother; plot signal power; make sure FIR size is correct, it was off by one

calibration/src/sphear_fir_design.m

@@ -1,10 +1,11 @@
+%% -*- mode: octave -*-
 %%
 %% *SpHEAR Project
 %%
 %% sphear_fir_design: design a matrix of FIR filters that equalizes the microphone
 %%
 %% usage:
-%%   [FIR] = sphear_fir_design(F, A2BF, n, Fs, Fl, plot_on, W, X, Y, Z)
+%%   [FIR] = sphear_fir_design(F, A2BF, n, maxboost, Fs, Fl, plot_on, W, X, Y, Z)
 %%
 %% F: array of filter section center frequencies
 %% A2BF: vector of A2B matrixes, one for each frequency
@@ -63,66 +64,122 @@ function [FIR] = sphear_fir_design(F, A2BF, n, maxboost, Fs, Fl, plot_on, W, X,
   %% create a frequency response profile in dB
   fr_profile = interp1([F(1), F(fl_idx), F(end)], [0, 0, 0], F);
 
+  %% boost high frequencies
+  [boost_l_B, boost_l_A] = boost(2.0, 8000, 15000, Fs);
+  [boost_h_B, boost_h_A] = boost(2.0, 12000, 10000, Fs);
+  [boost_l_h, boost_l_w] = freqz(boost_l_B, boost_l_A);
+  [boost_h_h, boost_h_w] = freqz(boost_h_B, boost_h_A);
+  %%fr_profile = interp1(boost_l_w/pi*(Fs/2), 20*log10(abs(boost_l_h)) .+ 20*log10(abs(boost_h_h)), F);
+
   dmax = size(A2BF(:,1,1))(1);
-  %% all b format components
+  %% calculate running averages of WXYZ filters and A2B matrices
   for b = 1:4
     %% all capsules
     for c = 1:4
-      %% all frequency bands
+      %% all frequency bands going down
+      for d = dmax:-1:1
+	if (F(d) > Fl)
+	  %% running average of A2B matrices
+	  if (d == dmax)
+	    upper_a2b(c, b, :) = A2BF(d, c, b);
+	  else
+	    if (d < dmax)
+	      upper_a2b(c, b, :) = (upper_a2b(c, b, :) + A2BF(d, c, b)) / 2;
+	    end
+	  end
+	  if (exist("W") == 1)
+	    %% running average of the filter values
+	    if (d == dmax)
+	      upper_w(c, b) = W(d);
+	      upper_x(c, b) = X(d);
+	      upper_y(c, b) = Y(d);
+	      upper_z(c, b) = Z(d);
+	    else
+	      if (d < dmax)
+		upper_w(c, b) = (upper_w(c, b) + W(d)) / 2;
+		upper_x(c, b) = (upper_x(c, b) + X(d)) / 2;
+		upper_y(c, b) = (upper_y(c, b) + Y(d)) / 2;
+		upper_z(c, b) = (upper_z(c, b) + Z(d)) / 2;
+	      end
+	    end
+	  end
+	end
+      end
+      %% all frequency bands going up
       for d = 1:dmax
-	%% create an array of points to design the filter
 	if (F(d) <= Fl)
-	  %% below the limit frequency we use the A2B coefficients
-	  vlast = A2BF(d, c, b);
-	  %% apply frequency response profile
-	  A2BF(d, c, b) = vlast * 10^(fr_profile(d)/20);
-	  dlast = d;
+	  %% running average of A2B matrices
+	  if (d == 1)
+	    lower_a2b(c, b, :) = A2BF(d, c, b);
+	  else
+	    if (d > 1)
+	      lower_a2b(c, b, :) = (lower_a2b(c, b, :) + A2BF(d, c, b)) / 2;
+	    end
+	  end
 	  if (exist("W") == 1)
+	    %% running average of the filter values
 	    if (d == 1)
-	      wlast = W(d);
-	      xlast = X(d);
-	      ylast = Y(d);
-	      zlast = Z(d);
+	      lower_w(c, b) = W(d);
+	      lower_x(c, b) = X(d);
+	      lower_y(c, b) = Y(d);
+	      lower_z(c, b) = Z(d);
 	    else
 	      if (d > 1)
-		wlast = (wlast + W(d)) / 2;
-		xlast = (xlast + X(d)) / 2;
-		ylast = (ylast + Y(d)) / 2;
-		zlast = (zlast + Z(d)) / 2;
+		lower_w(c, b) = (lower_w(c, b) + W(d)) / 2;
+		lower_x(c, b) = (lower_x(c, b) + X(d)) / 2;
+		lower_y(c, b) = (lower_y(c, b) + Y(d)) / 2;
+		lower_z(c, b) = (lower_z(c, b) + Z(d)) / 2;
 	      end
 	    end
 	  end
+	end
+      end
+    end
+  end
+
+  %% all b format components
+  for b = 1:4
+    %% all capsules
+    for c = 1:4
+      %% all frequency bands
+      for d = 1:dmax
+	%% create an array of points to design the filter
+	if (F(d) <= Fl)
+	  dlast = d;
+	  %% apply frequency response profile
+	  A2BF(d, c, b) = A2BF(d, c, b) * 10^(fr_profile(d)/20);
 	else
 	  %% above the limit frequency we insert the W, X, Y and Z
 	  %% data if they exist, otherwise keep the last valid value
 	  %% from the A2B coefficients (ie: flat above the transition
 	  %% frequency
+	  limit_a2b = ((lower_a2b(c, b, :) + lower_a2b(c, b, :)) / 2);
 	  if (b == 1)
 	    if (exist("W") == 1)
-	      A2BF(d, c, b) = vlast * (W(d)/wlast) * 10^(fr_profile(d)/20);
+	      A2BF(d, c, b) = limit_a2b * (W(d)/((lower_w(c, b) + upper_w(c, b)) / 2)) * 10^(fr_profile(d)/20);
 	    else
-	      A2BF(d, c, b) = vlast;
+	      A2BF(d, c, b) = limit_a2b;
 	    end
 	  end
 	  if (b == 2)
 	    if (exist("X") == 1)
-	      A2BF(d, c, b) = vlast * (X(d)/xlast) * 10^(fr_profile(d)/20);
+	      A2BF(d, c, b) = limit_a2b * (X(d)/((lower_x(c, b) + upper_x(c, b)) / 2)) * 10^(fr_profile(d)/20);
 	    else
-	      A2BF(d, c, b) = vlast;
+	      A2BF(d, c, b) = limit_a2b;
 	    end
 	  end
 	  if (b == 3)
 	    if (exist("Y") == 1)
-	      A2BF(d, c, b) = vlast * (Y(d)/ylast) * 10^(fr_profile(d)/20);
+	      A2BF(d, c, b) = limit_a2b * (Y(d)/((lower_y(c, b) + upper_y(c, b)) / 2)) * 10^(fr_profile(d)/20);
 	    else
-	      A2BF(d, c, b) = vlast;
+	      A2BF(d, c, b) = limit_a2b;
 	    end
 	  end
 	  if (b == 4)
 	    if (exist("Z") == 1)
-	      A2BF(d, c, b) = vlast * (Z(d)/zlast) * 10^(fr_profile(d)/20);
+	      A2BF(d, c, b) = limit_a2b * (Z(d)/((lower_z(c, b) + upper_z(c, b)) / 2)) * 10^(fr_profile(d)/20);
 	    else
-	      A2BF(d, c, b) = vlast;
+	      A2BF(d, c, b) = limit_a2b;
 	    end
 	  end
 	end
@@ -164,8 +221,8 @@ function [FIR] = sphear_fir_design(F, A2BF, n, maxboost, Fs, Fl, plot_on, W, X,
 
       %% we fit the spline taking into account one additional point
       %% on the bottom edge so that it will fit better with the existing
-      %% low frequency points, now we need to find the exact transition
-      %% frequency in the interpolated array
+      %% low frequency points, now we find the exact transition frequency
+      %% in the interpolated array
       for idx = 1:size(fh(:))
 	if (fh(idx) <= F(dlast)/(Fs/2))
 	  istart = idx;
@@ -176,22 +233,21 @@ function [FIR] = sphear_fir_design(F, A2BF, n, maxboost, Fs, Fl, plot_on, W, X,
       fl = [0 F(1:dlast-1)']' ./ (Fs/2);
       rl = [A2BF(1,c,b) A2BF(1:dlast-1,c,b)']';
 
+%%      fh = [F(dlast:end)' Fs/2]' ./ (Fs/2);
+%%      rh = [A2BF(dlast:end,c,b)' A2BF(end,c,b)]';
+
       %% concatenate low and high frequency arrays
       fi = cat(1, fl, fh(istart:end));
       ri = cat(1, rl, rh(istart:end));
 
-      printf("low frequency boost %f, high frequency boost %f\n", 20*log10(ri(1)), 20*log10(ri(end)));
       %% force a max boost at low and high frequencies
       if (maxboost != 0)
-	if (abs(ri(1)) > 10^(maxboost/20))
-	  printf("limit low freq boost from %f to %f dB (%f)\n", 20*log10(ri(1)), maxboost, ri(1));
-	  ri(1) = sign(ri(1)) * 10^(maxboost/20);
-	  printf("new limit %f\n", ri(1));
-	end
-	if (abs(ri(end)) > 10^(maxboost/20))
-	  printf("limit high freq boost from %f dB to %f dB (%f)\n", 20*log10(ri(end)), maxboost, ri(end));
-	  ri(end) = sign(ri(end)) * 10^(maxboost/20);
-	  printf("new limit %f\n", ri(end));
+	printf("b-format component:%d capsule:%d \n", b, c);
+	for idx = 1:size(ri(:))
+	  if (abs(ri(idx)) > 10^(maxboost/20))
+	    printf("freq=%.2f: limit boost from %.2f to %.2f dB\n", fi(idx)*Fs/2, 10*log10(ri(idx)), maxboost);
+	    ri(idx) = sign(ri(idx)) * 10^(maxboost/20);
+	  end
 	end
       end
       
@@ -224,18 +280,18 @@ function [FIR] = sphear_fir_design(F, A2BF, n, maxboost, Fs, Fl, plot_on, W, X,
 	%% title(sprintf("B format FIR filters (transition frequency: %d Hz)", Fl));
 	hold on;
 	%% plot real measured points
-	%% ph = semilogx(F, 20*log10(abs(A2BF(:,c,b)))', popts, "linewidth", 1.5);
+	%% ph = semilogx(F, 10*log10(abs(A2BF(:,c,b)))', popts, "linewidth", 1.5);
 	%% plot extrapolated filter
 	if (b == 1)
 	  %% W
-	  wh = semilogx(fi*(Fs/2), 20*log10(abs(ri)), wopts, "linewidth", 1.5);
+	  wh = semilogx(fi*(Fs/2), 10*log10(abs(ri)), wopts, "linewidth", 1.5);
 	else
 	  if (b == 4)
 	    %% Z
-	    zh = semilogx(fi*(Fs/2), 20*log10(abs(ri)), zopts, "linewidth", 1.5);
+	    zh = semilogx(fi*(Fs/2), 10*log10(abs(ri)), zopts, "linewidth", 1.5);
 	  else
 	    %% XY
-	    xyh = semilogx(fi*(Fs/2), 20*log10(abs(ri)), xyopts, "linewidth", 1.5);
+	    xyh = semilogx(fi*(Fs/2), 10*log10(abs(ri)), xyopts, "linewidth", 1.5);
 	  end
 	end
       end
@@ -245,7 +301,8 @@ function [FIR] = sphear_fir_design(F, A2BF, n, maxboost, Fs, Fl, plot_on, W, X,
       %% "Two alternative minimum-phase filters tested perceptually"
       %% Robert Mores and Ralf Hendrych, AES Convertion paper 9554
       %%
-      IR = fir2(n, fi, ri)';
+      IR = fir2(n - 1, fi, ri)';
+      %%
       %% number of samples in the sin transition
       %%
       %% 8 seems to be the sweet spot for a FIR filter size of 256,
@@ -253,22 +310,26 @@ function [FIR] = sphear_fir_design(F, A2BF, n, maxboost, Fs, Fl, plot_on, W, X,
       %% the flattness of the frequency response
       %%
       M = size(IR)(1);
-      N = (M-1)/32;
+      N = M/32;
       R = 0.5 + (0.5 * sin(pi * ((1:(2*N)) - N) / (2 * N)));
-      ramp = cat(1, zeros(M/2 - N + 1, 1), R', ones(M/2 - N, 1));
+      ramp = cat(1, zeros(M/2 - N, 1), R', ones(M/2 - N, 1));
       IRMP = IR .* ramp * 2;
-
+      %%
       %% now we trim the filter and return one that is 1/2 the size requested
+      %% because it is now almost minimum phase and the first half is zeroed
+      %%
       %% we also window the end of the impulse response... (unnecessary?)
       %% T = 4;
       %% tapper = cat(1, ones(M/2 - 2*T + 1, 1), (1.0 - (0.5 + (0.5 * sin(pi * ((1:(2*T)) - T) / (2 * T)))))');
-      FIR(b, c, :) = IRMP((M-1)/2 - N:(M-1)/2 - N + (M-1)/2);
+      %%
+      FIR(b, c, :) = IRMP(M/2 - N:M/2 - N + M/2 - 1);
     end
   end
 
   %% write a tetraproc compatible wav file
   fuma = [sqrt(2)/2, 1, 1, 1];
   for b = 1:4
+    %% each column is one channel, W/X/Y/Z, all capsules concatenated
     TETRA(:, b) = fuma(b) * cat(1,
 				reshape(FIR(b, 1, :), 1, [])',
 				reshape(FIR(b, 2, :), 1, [])',