calculate filters for X Y and Z separately, no need to average all components at once

sphear_fir_design.m

@@ -4,21 +4,23 @@
 %% sphear_fir_design: design a matrix of FIR filters that equalizes the microphone
 %%
 %% usage:
-%%   [FIR] = sphear_fir_design(F, A2BF, n, Fs, Fl, W, XYZ)
+%%   [FIR] = sphear_fir_design(F, A2BF, n, Fs, Fl, W, X, Y, Z)
 %%
 %% 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
-%% Fl: frequency limit, below it use A2B matrixes, above it use WXYZ shapes
+%% Fl: frequency limit, below it use A2B matrixes, above it use W, X, Y, Z shapes
 %%
-%% for the first run do not supply W/XYZ, render the filters and then run
-%% sphear_fir_extract to extract the W and XYZ filter shapes, finally run
+%% 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
 %%
 %% W: extracted shape of W filter (used above Fl)
-%% XYZ: extracted shape of X, Y and Z filters (used above Fl)
+%% X: extracted shape of X filter (used above Fl)
+%% Y: extracted shape of Y filter (used above Fl)
+%% Z: extracted shape of Z filter (used above Fl)
 %%
 %% Copyright 2016, Fernando Lopez-Lezcano, All Rights Reserved
 %% nando@ccrma.stanford.edu
@@ -36,7 +38,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, W, XYZ)
+function [FIR] = sphear_fir_design(F, A2BF, n, maxboost, Fs, Fl, W, X, Y, Z)
 
   fig = figure;
 
@@ -52,12 +54,12 @@ function [FIR] = sphear_fir_design(F, A2BF, n, maxboost, Fs, Fl, W, XYZ)
 	  dlast = d;
 	  if (exist("W") == 1)
 	    wlast = W(d);
-	  end
-	  if (exist("XYZ") == 1)
-	    xyzlast = XYZ(d);
+	    xlast = X(d);
+	    ylast = Y(d);
+	    zlast = Z(d);
 	  end
 	else
-	  %% above the limit frequency we interpolate the W and XYZ
+	  %% above the limit frequency we interpolate the W, X, Y and Z
 	  %% filters if they exist, otherwise keep the last valid value
 	  if (b == 1)
 	    if (exist("W") == 1)
@@ -65,9 +67,24 @@ function [FIR] = sphear_fir_design(F, A2BF, n, maxboost, Fs, Fl, W, XYZ)
 	    else
 	      A2BF(d, c, b) = vlast;
 	    end
-	  else
-	    if (exist("XYZ") == 1)
-	      A2BF(d, c, b) = vlast * (XYZ(d)/xyzlast);
+	  end
+	  if (b == 2)
+	    if (exist("X") == 1)
+	      A2BF(d, c, b) = vlast * (X(d)/xlast);
+	    else
+	      A2BF(d, c, b) = vlast;
+	    end
+	  end
+	  if (b == 3)
+	    if (exist("Y") == 1)
+	      A2BF(d, c, b) = vlast * (Y(d)/ylast);
+	    else
+	      A2BF(d, c, b) = vlast;
+	    end
+	  end
+	  if (b == 4)
+	    if (exist("Z") == 1)
+	      A2BF(d, c, b) = vlast * (Z(d)/zlast);
 	    else
 	      A2BF(d, c, b) = vlast;
 	    end
@@ -75,7 +92,7 @@ function [FIR] = sphear_fir_design(F, A2BF, n, maxboost, Fs, Fl, W, XYZ)
 	end
       end
 
-      %% fit a spline to the W and XYZ segments and extrapolate
+      %% fit a spline to the W, X, Y and Z segments and extrapolate
       %% highest frequency if it does not exist in the measurements
       %%
       %% the 200Hz spacing of the points is arbitrary, we should

sphear_fir_extract.m

@@ -22,7 +22,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 [W XYZ] = sphear_fir_extract(B, all, onaxis);
+function [W X Y Z] = sphear_fir_extract(B, all, onaxis);
 
   if (all == 1)
     m = (1:16);
@@ -38,7 +38,11 @@ function [W XYZ] = sphear_fir_extract(B, all, onaxis);
   W = 1.0 ./ reshape(mean(B(m,1,:)), size(B(:,:,:))(3), 1);
 
   %% extract average response shape for XYZ
-  XYZ = 1.0 ./ reshape(mean(mean(B(m, [2, 3, 4], :))), size(B(:,:,:))(3), 1);
+  % XYZ = 1.0 ./ reshape(mean(mean(B(m, [2, 3, 4], :))), size(B(:,:,:))(3), 1);
+
+  X = 1.0 ./ reshape(mean(B(m,2,:)), size(B(:,:,:))(3), 1);
+  Y = 1.0 ./ reshape(mean(B(m,3,:)), size(B(:,:,:))(3), 1);
+  Z = 1.0 ./ reshape(mean(B(m,4,:)), size(B(:,:,:))(3), 1);
 
 end