remove temporary TFIR return parameter, some changes to parameters to function

calibration/src/sphear_design_fir_filters.m

@@ -2,31 +2,28 @@
 %%
 %% *SpHEAR Project
 %%
-%% sphear_design_fir_filters: design FIR filters based on W, X, Y and Z shapes
+%% sphear_design_fir_filters: design FIR filters based on B format spectrum shapes
 %%
 %% usage:
-%%   [FIR] = sphear_design_fir_filters(F, n, maxboost, Fs, plot_on,
-%%                                     shelf_gain, shelf_freq, W, X, Y, Z)
+%%   [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)
 %%
 %% F: array of filter section center frequencies
 %% n: order of the FIR filter being created
-%% maxboost: max high and low frequency boost in dB of the final filter
 %% Fs: sampling rate
 %% plot_on: plot frequency response of components, default is "epsc"
 %%          can be one of 0, "eps", "epsc" (color), "png", "pngc" (color)
 %%
-%% shelf_gain: optional high pass shelf filter, gain in dB (0.0 dB turns it off)
-%% shelf_freq: corner frequency in Hertz
+%% 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
 %%
-%% 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
+%% LF_max_boost_db: maximum boost of the second order compensation filters
+%% DC_gain_offset_db: adjust gain at DC by this offset
 %%
-%% W: extracted shape of W filter (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)
-%%    Z is optional, if not supplied X and Y are averaged to obtain Z
+%% FIR: array of minimum phase FIR filters
+%% IR: calculated FIR filters
+%% MPS: output of mps without windowing (minimum phase)
 %%
 %% Copyright 2016-2017, Fernando Lopez-Lezcano, All Rights Reserved
 %% nando@ccrma.stanford.edu
@@ -44,29 +41,30 @@
 %%  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 TFIR IR] = sphear_design_fir_filters(F, n, maxboost, SR, W, X, Y, Z, U, V, LF_max_boost, DC_gain_offset)
+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)
   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);
+  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")
-    %% define gain offset at DC for UV filters
-    DC_gain_offset = db2mag(-30);
+  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 because of the minimal
-  %% phase calculation
+  %% 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)]');
-  if (exist("U"))
+  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)]');
   end
@@ -78,14 +76,14 @@ function [FIR TFIR IR] = sphear_design_fir_filters(F, n, maxboost, SR, W, X, Y,
       %%
       %% https://ccrma.stanford.edu/~jos/filters/Conversion_Minimum_Phase.html
       %%
-      FIR(b,:) = real(ifft(mps(fft(IR(b,:)))));
+      MPS(b,:) = real(ifft(mps(fft(IR(b,:)))));
       %% create a raised cosine window for the end of the IR
-      m = n/2;
+      m = n;
       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;
+      FIR(b,:) = (MPS(b,:)')(1:m) .* end_window;
     else
-      %% try to get the filter to be close to minimum phase
+      %% 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
@@ -100,7 +98,7 @@ function [FIR TFIR IR] = sphear_design_fir_filters(F, n, maxboost, SR, W, X, Y,
       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;
+      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
@@ -109,9 +107,7 @@ function [FIR TFIR IR] = sphear_design_fir_filters(F, n, maxboost, SR, W, X, Y,
       %% 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,:);
+      FIR(b,:) = MPS(M/2 - N:M/2 - N + M/2 - 1);
     end
   end
 end