sphear_design_fir_filters: fix n-1 mistake when making filters, fix linear phase code to match, remove old code

calibration/src/sphear_design_fir_filters.m

@@ -5,21 +5,21 @@
 %% sphear_design_fir_filters: design FIR filters based on B format spectrum shapes
 %%
 %% usage:
-%%   [FIR IR MPS] = sphear_design_fir_filters(F, n, Fs, plot_on, W, X, Y, Z, U, V,
-%%                                            LF_max_boost_db, DC_gain_offset_db)
+%%   [FIR IR MPS] = sphear_design_fir_filters(F, n, SR, W, X, Y, Z, U, V,
+%%                                            UV_lf_freq, UV_lf_slope, DC_gain_offset_db)
 %%
 %% F: array of filter section center frequencies
 %% n: order of the FIR filter being created
-%% Fs: sampling rate
-%% plot_on: plot frequency response of components, default is "epsc"
-%%          can be one of 0, "eps", "epsc" (color), "png", "pngc" (color)
+%% SR: sampling rate
 %%
 %% W, X, Y, Z, U, V: extracted shape of B format component filter (used above Fl)
 %%    Z: optional, if not supplied X and Y are averaged to obtain Z
 %%    UV: optional, if not supplied only first order shapes are calculated
+%%        can be arrays and then an array of filters will be calculated and returned
 %%
-%% LF_max_boost_db: maximum boost of the second order compensation filters
-%% DC_gain_offset_db: adjust gain at DC by this offset
+%% UV_lf_freq: corner frequency of low pass filter that limits low frequency gain of UV
+%% UV_lf_slop: slope in dB/oct of the low pass filter
+%% DC_gain_offset_db: DC filter gain in dB
 %%
 %% FIR: array of minimum phase FIR filters
 %% IR: calculated FIR filters
@@ -41,74 +41,66 @@
 %%  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 IR MPS] = sphear_design_fir_filters(F, n, SR, W, X, Y, Z, U, V, LF_max_boost_db, DC_gain_offset_db)
+%% NOTES: using fir2 for now, firls calculates non-smooth filters with worse
+%% frequency response, remez does not converge in my tests
+
+function [FIR IR MPS] = sphear_design_fir_filters(F, n, SR, W, X, Y, Z, U, V, UV_lf_freq, UV_lf_slope, DC_gain_offset_db)
+  %% arbitrary oversampling for FIR length
   if ! exist("Z")
     %% if Z was not supplied, just average X and Y data
     Z = (X + Y) / 2;
   end
-  if ! exist("LF_max_boost_db")
-    %% by default no maximum UV boost at low frequencies
-    LF_max_boost_db = 100;
-  end
   if ! exist("DC_gain_offset_db")
     %% by default no gain offset at DC for UV filters
     DC_gain_offset_db = 0;
   end
-  LF_max_boost = db2mag(LF_max_boost_db);
   DC_gain_offset = db2mag(DC_gain_offset_db);
   %% 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
   %% final filters will become 1/2 the size when folded for minimum phase
-  IR(1,:) = fir2((n * 2) - 1, freqs, [W(1), W', W(end)]');
-  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)]');
+  IR(1,:) = fir2(n, freqs, [W(1), W', W(end)]');
+  IR(2,:) = fir2(n, freqs, [X(1), X', X(end)]');
+  IR(3,:) = fir2(n, freqs, [Y(1), Y', Y(end)]');
+  IR(4,:) = fir2(n, freqs, [Z(1), Z', Z(end)]');
   if (exist("U") && exist("V"))
-    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)]');
+    %% make a n-dB/oct low pass filter shape for UV components
+    %% this limits the high pass filter gain (and noise) at low frequencies
+    for f = (1:size(F))
+      if (F(f) > UV_lf_freq)
+	UV_lpf(f,:) = 1.0;
+      else
+	UV_lpf(f,:) = db2mag((F(f) - UV_lf_freq)/UV_lf_freq*UV_lf_slope);
+      end
+    end
+    UV_lpf = [DC_gain_offset, UV_lpf', 1.0]';
+    %% calculate UV filters
+    if (size(U, 2) > 1)
+      %% we are using one U/V filter for each capsule
+      %% add all filters at end of array (16 total)
+      for c = 1:size(U, 2)
+	IR(4+c,:) = fir2(n, freqs, [U(1,c), U(:,c)', U(end,c)]' .* UV_lpf);
+	IR(4+size(U,2)+c,:) = fir2(n, freqs, [V(1,c), V(:,c)', V(end,c)]' .* UV_lpf);
+      end
+    else
+      IR(5,:) = fir2(n, freqs, [U(1), U', U(end)]' .* UV_lpf);
+      IR(6,:) = fir2(n, freqs, [V(1), V', V(end)]' .* UV_lpf);
+    end
   end
 
   for b = (1:size(IR, 1))
-    if (1)
-      %%
-      %% convert FIR filter to minimum phase filter
-      %%
-      %% https://ccrma.stanford.edu/~jos/filters/Conversion_Minimum_Phase.html
-      %%
-      MPS(b,:) = real(ifft(mps(fft(IR(b,:)))));
-      %% create a raised cosine window for the end of the IR
-      m = n;
-      end_window = tukeywin(2 * m, 0.25)(m+1:end);
-      %% truncate to original size and window end
-      FIR(b,:) = (MPS(b,:)')(1:m) .* end_window;
-    else
-      %% try to get the filter to be close to minimum phase (obsolete code)
-      %%
-      %% "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));
-      MPS = 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,:) = MPS(M/2 - N:M/2 - N + M/2 - 1);
-    end
+    %%
+    %% convert FIR filter to minimum phase filter
+    %%
+    %% https://ccrma.stanford.edu/~jos/filters/Conversion_Minimum_Phase.html
+    %%
+    MPS(b,:) = real(ifft(mps(fft(IR(b,:)))));
+    %%
+    %% create a raised cosine window for the end of the IR
+    %% to make sure IR decays to zero (last 8 samples of a 256 sample IR)
+    end_window = tukeywin(2 * n, 0.025)(n+1:end);
+    %% truncate to original size and window end
+    FIR(b,:) = (MPS(b,:)')(1:n) .* end_window;
   end
 end