add maximum amount of boost to limit filter ends after fitting the spline, this still does not work out in low frequencies, perhaps we need to fit the spline only in the upper frequency range and just leave the A2B matrix alone in the low frequency range

sphear_fir_design.m

@@ -9,6 +9,7 @@
 %% F: array of filter section center frequencies
 %% A2BF: vector of A2B matrixes, one for each frequency
 %% n: order of the FIR filter being created
+%% maxboost: max high and low frequency boost in dB of the final filter
 %% Fs: sampling rate
 %%
 %% for the first run do not supply Fl/W/XYZ, render the filters and from inspection
@@ -35,13 +36,10 @@
 %%  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, Fs, Fl, W, XYZ)
+function [FIR] = sphear_fir_design(F, A2BF, n, maxboost, Fs, Fl, W, XYZ)
 
   FIR = zeros(4, 4, (n + 1));
 
-  %% frequency range
-  f = [0 F' Fs/2]' / (Fs/2);
-
   %% all b format components
   for b = 1:4
     %% all capsules
@@ -59,7 +57,8 @@ function [FIR] = sphear_fir_design(F, A2BF, n, Fs, Fl, W, XYZ)
 	    xyzlast = XYZ(d);
 	  end
 	else
-	  %% above the limit frequency we interpolate the W and XYZ shapes
+	  %% above the limit frequency we interpolate the W and XYZ
+	  %% filters if they exist, otherwise keep the last valid value
 	  if (b == 1)
 	    if (exist("W") == 1)
 	      A2BF(d, c, b) = vlast * (W(d)/wlast);
@@ -75,27 +74,41 @@ function [FIR] = sphear_fir_design(F, A2BF, n, Fs, Fl, W, XYZ)
 	  end
 	end
       end
-      
-      %% interpolate last point at nyquist from the previous two points
-      r = [ A2BF(1,c,b) A2BF(:,c,b)' (A2BF(end,c,b) + (A2BF(end,c,b) - A2BF(end-1,c,b)))]';
 
-      %% get a spline fitted over the points to smooth the response
-      %% endpoints have null first derivative
+      %% fit a spline over the existing points and extrapolate
+      %% for low and high ends of spectrum
       fi = (0:200/(Fs/2):1)';
-      ri = spline(f, [0 r' 0], fi);
+      ri = spline(F/(Fs/2), A2BF(:,c,b)', fi);
 
+      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));
+	end
+      end
+      
+      %% real measured points
       if (b == 1)
-	semilogx(f, r,'*g');
+	semilogx(F/(Fs/2), A2BF(:,c,b)','*g');
       else
-	semilogx(f, r,'*b');
+	semilogx(F/(Fs/2), A2BF(:,c,b)','*b');
       end
       hold on;
+      %% extrapolated spline
       if (b == 1)
 	semilogx(fi,ri,'g');
       else
 	semilogx(fi,ri,'r');
       end
-      
+
       %% design the filter
       FIR(b, c, :) = fir2(n, fi, ri);
     end