calibration: add optional high pass shelf filter to the design function, this can be used to enhance the high frequency response of the microphone

calibration/src/sphear_fir_design.m

@@ -5,7 +5,8 @@
 %% sphear_fir_design: design a matrix of FIR filters that equalizes the microphone
 %%
 %% usage:
-%%   [FIR] = sphear_fir_design(F, A2BF, n, maxboost, Fs, Fl, plot_on, W, X, Y, Z)
+%%   [FIR] = sphear_fir_design(F, A2BF, n, maxboost, Fs, Fl, plot_on,
+%%                             shelf_gain, shelf_freq, transition_octaves, W, X, Y, Z)
 %%
 %% F: array of filter section center frequencies
 %% A2BF: vector of A2B matrixes, one for each frequency
@@ -16,6 +17,10 @@
 %% plot_on: plot frequency response of components, default is "epsc"
 %%          can be one of 0, "eps", "epsc" (color), "png", "pngc" (color)
 %%
+%% shelf_gain: optional high pass shelf filter, gain in dB (0.0 dB turns it off)
+%% shelf_freq: corner frequency in Hertz
+%% transition_octaves: transition width in octaves below the corner frequency
+%%
 %% for the first run do not supply W/X/Y/Z, render the filters and then run
 %% sphear_fir_extract to extract the W, X, Y and Z filter shapes, finally run
 %% this again with the extra parameters below
@@ -42,7 +47,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 [FIR] = sphear_fir_design(F, A2BF, n, maxboost, Fs, Fl, plot_on, W, X, Y, Z)
+function [FIR] = sphear_fir_design(F, A2BF, n, maxboost, Fs, Fl, plot_on, shelf_gain, shelf_freq, transition_octaves, W, X, Y, Z)
 
   if (plot_on != 0)
     fig = figure;
@@ -61,15 +66,17 @@ function [FIR] = sphear_fir_design(F, A2BF, n, maxboost, Fs, Fl, plot_on, W, X,
     end
   end
   printf("transition band is %d at %f Hz\n", fl_idx, F(fl_idx));
-  %% 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);
+  %% create optional high pass shelf filter (in dB)
+  if (exist("shelf_gain") && (shelf_gain != 0.0))
+    shelf_corner = shelf_freq / (2^transition_octaves);
+    %% transition
+    shelf = (tukeywin(32, 0.9)(1:end/2) * shelf_gain);
+    %% interpolate over existing frequency array
+    fr_profile = interp1([F(1), linspace(shelf_corner, shelf_freq, 16), F(end)], [0.0, shelf', shelf_gain], F);
+  else
+    fr_profile = interp1([F(1), F(fl_idx), F(end)], [0, 0, 0], F);
+  end
 
   dmax = size(A2BF(:,1,1))(1);
   %% calculate running averages of WXYZ filters and A2B matrices
@@ -147,7 +154,7 @@ function [FIR] = sphear_fir_design(F, A2BF, n, maxboost, Fs, Fl, plot_on, W, X,
 	if (F(d) <= Fl)
 	  dlast = d;
 	  %% apply frequency response profile
-	  A2BF(d, c, b) = A2BF(d, c, b) * 10^(fr_profile(d)/20);
+	  A2BF(d, c, b) = A2BF(d, c, b) * 10^(fr_profile(d)/10);
 	else
 	  %% above the limit frequency we insert the W, X, Y and Z
 	  %% data if they exist, otherwise keep the last valid value
@@ -156,28 +163,28 @@ function [FIR] = sphear_fir_design(F, A2BF, n, maxboost, Fs, Fl, plot_on, W, X,
 	  limit_a2b = ((lower_a2b(c, b, :) + lower_a2b(c, b, :)) / 2);
 	  if (b == 1)
 	    if (exist("W") == 1)
-	      A2BF(d, c, b) = limit_a2b * (W(d)/((lower_w(c, b) + upper_w(c, b)) / 2)) * 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)/10);
 	    else
 	      A2BF(d, c, b) = limit_a2b;
 	    end
 	  end
 	  if (b == 2)
 	    if (exist("X") == 1)
-	      A2BF(d, c, b) = limit_a2b * (X(d)/((lower_x(c, b) + upper_x(c, b)) / 2)) * 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)/10);
 	    else
 	      A2BF(d, c, b) = limit_a2b;
 	    end
 	  end
 	  if (b == 3)
 	    if (exist("Y") == 1)
-	      A2BF(d, c, b) = limit_a2b * (Y(d)/((lower_y(c, b) + upper_y(c, b)) / 2)) * 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)/10);
 	    else
 	      A2BF(d, c, b) = limit_a2b;
 	    end
 	  end
 	  if (b == 4)
 	    if (exist("Z") == 1)
-	      A2BF(d, c, b) = limit_a2b * (Z(d)/((lower_z(c, b) + upper_z(c, b)) / 2)) * 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)/10);
 	    else
 	      A2BF(d, c, b) = limit_a2b;
 	    end
@@ -241,9 +248,9 @@ function [FIR] = sphear_fir_design(F, A2BF, n, maxboost, Fs, Fl, plot_on, W, X,
       if (maxboost != 0)
 	printf("b-format component:%d capsule:%d \n", b, c);
 	for idx = 1:size(ri(:))
-	  if (abs(ri(idx)) > 10^(maxboost/20))
+	  if (abs(ri(idx)) > 10^(maxboost/10))
 	    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);
+	    ri(idx) = sign(ri(idx)) * 10^(maxboost/10);
 	  end
 	end
       end