sphear_fir_design: change completely to use supplied W and XYZ filter shapes for the high frequency region

sphear_fir_design.m

@@ -4,13 +4,21 @@
 %% sphear_fir_design: design a matrix of FIR filters that equalizes the microphone
 %%
 %% usage:
-%%   [FIR] = sphear_fir_design(F, A2BF, n, Fs)
+%%   [FIR] = sphear_fir_design(F, A2BF, n, Fs, Fl, W, XYZ)
 %%
 %% F: array of filter section center frequencies
-%% A2BF: vector of A2B matrixes
+%% A2BF: vector of A2B matrixes, one for each frequency
 %% n: order of the FIR filter being created
 %% Fs: sampling rate
 %%
+%% for the first run do not supply Fl/W/XYZ, render the filters and from inspection
+%% determine the limit frequency, then run sphear_fir_extract to extract the
+%% W and XYZ filter shapes and run this again with the extra parameters below
+%%
+%% Fl: frequency limit, below it use A2B matrixes, above it use WXYZ shapes
+%% W: extracted shape of W filter (used above Fl)
+%% XYZ: extracted shape of X, Y and Z filters (used above Fl)
+%%
 %% Copyright 2016, Fernando Lopez-Lezcano, All Rights Reserved
 %% nando@ccrma.stanford.edu
 %%
@@ -27,7 +35,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, Fs)
+function [FIR] = sphear_fir_design(F, A2BF, n, Fs, Fl, W, XYZ)
 
   FIR = zeros(4, 4, (n + 1));
 
@@ -38,7 +46,37 @@ function [FIR] = sphear_fir_design(F, A2BF, n, Fs)
   for b = 1:4
     %% all capsules
     for c = 1:4
+
+      for d = 1:size(A2BF(:,1,1))(1)
+	if (F(d) <= Fl)
+	  %% below the limit frequency we use the A2B coefficients
+	  vlast = A2BF(d, c, b);
+	  if (exist("W") == 1)
+	    wlast = W(d);
+	  end
+	  if (exist("XYZ") == 1)
+	    xyzlast = XYZ(d);
+	  end
+	else
+	  %% above the limit frequency we interpolate the W and XYZ shapes
+	  if (b == 1)
+	    if (exist("W") == 1)
+	      A2BF(d, c, b) = vlast * (W(d)/wlast);
+	    else
+	      A2BF(d, c, b) = vlast;
+	    end
+	  else
+	    if (exist("XYZ") == 1)
+	      A2BF(d, c, b) = vlast * (XYZ(d)/xyzlast);
+	    else
+	      A2BF(d, c, b) = vlast;
+	    end
+	  end
+	end
+      end
+      %% full response
       r = [ A2BF(1,c,b) A2BF(:,c,b)' A2BF(end,c,b)]';
+      %% design the filter
       FIR(b, c, :) = fir2(n, f, r);
     end
   end