transform FIR filters into proper minimum phase filters, add support for second order horizontal components and preliminary low frequency gain limiting of the second order filters

calibration/src/sphear_design_fir_filters.m

@@ -44,11 +44,19 @@
 %%  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_design_fir_filters(F, n, maxboost, SR, W, X, Y, Z)
+function [FIR TFIR IR] = sphear_design_fir_filters(F, n, maxboost, SR, W, X, Y, Z, U, V, LF_max_boost, DC_gain_offset)
   if ! exist("Z")
     %% if Z was not supplied, just average X and Y data
     Z = (X + Y) / 2;
   end
+  if ! exist("LF_max_boost")
+    %% define default maximum UV boost at low frequencies
+    LF_max_boost = db2mag(8);
+  end
+  if ! exist("DC_gain_offset")
+    %% define gain offset at DC for UV filters
+    DC_gain_offset = db2mag(-30);
+  end
   %% make a 0 to 1 vector of relative frequencies
   freqs = [0, (F/(SR/2))', 1]';
   %% design the filters, we use twice the requested size as the
@@ -58,37 +66,52 @@ function [FIR] = sphear_design_fir_filters(F, n, maxboost, SR, W, X, Y, Z)
   IR(2,:) = fir2((n * 2) - 1, freqs, [X(1), X', X(end)]');
   IR(3,:) = fir2((n * 2) - 1, freqs, [Y(1), Y', Y(end)]');
   IR(4,:) = fir2((n * 2) - 1, freqs, [Z(1), Z', Z(end)]');
+  if (exist("U"))
+    IR(5,:) = fir2((n * 2) - 1, freqs, [U(1) * DC_gain_offset, min(LF_max_boost, U)', U(end)]');
+    IR(6,:) = fir2((n * 2) - 1, freqs, [V(1) * DC_gain_offset, min(LF_max_boost, V)', V(end)]');
+  end
 
-  for b = (1:4)
+  for b = (1:size(IR, 1))
     if (1)
-    %% try to get the filter to be close to minimum phase
-    %%
-    %% "Two alternative minimum-phase filters tested perceptually"
-    %% Robert Mores and Ralf Hendrych, AES Convention paper 9554
-    %%
-    %% number of samples in the sin transition
-    %%
-    %% 8 seems to be the sweet spot for a FIR filter size of 256,
-    %% 16 leads to bigger phase shift, 4 starts to impact negatively
-    %% the flattness of the frequency response
-    %%
-    M = size(IR)(2);
-    N = M/32;
-    R = 0.5 + (0.5 * sin(pi * ((1:(2*N)) - N) / (2 * N)));
-    ramp = cat(1, zeros(M/2 - N, 1), R', ones(M/2 - N, 1));
-    IRMP = IR(b,:)' .* 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,:) = IRMP(M/2 - N:M/2 - N + M/2 - 1);
+      %%
+      %% convert FIR filter to minimum phase filter
+      %%
+      %% https://ccrma.stanford.edu/~jos/filters/Conversion_Minimum_Phase.html
+      %%
+      FIR(b,:) = real(ifft(mps(fft(IR(b,:)))));
+      %% create a raised cosine window for the end of the IR
+      m = n/2;
+      end_window = tukeywin(2 * m, 0.25)(m+1:end);
+      %% truncate to original size and window end
+      TFIR(b,:) = (FIR(b,:)')(1:m) .* end_window;
     else
-      size(IR(b,:))
-      FIR(b,:) = mps(fft(IR(b,:)));
+      %% try to get the filter to be close to minimum phase
+      %%
+      %% "Two alternative minimum-phase filters tested perceptually"
+      %% Robert Mores and Ralf Hendrych, AES Convention paper 9554
+      %%
+      %% number of samples in the sin transition
+      %%
+      %% 8 seems to be the sweet spot for a FIR filter size of 256,
+      %% 16 leads to bigger phase shift, 4 starts to impact negatively
+      %% the flattness of the frequency response
+      %%
+      M = size(IR)(2);
+      N = M/32;
+      R = 0.5 + (0.5 * sin(pi * ((1:(2*N)) - N) / (2 * N)));
+      ramp = cat(1, zeros(M/2 - N, 1), R', ones(M/2 - N, 1));
+      IRMP = IR(b,:)' .* 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,:) = IRMP(M/2 - N:M/2 - N + M/2 - 1);
+      %% XXX for testing
+      TFIR(b,:) = FIR(b,:);
     end
   end
 end