# HG changeset patch
# User dicklyon@google.com
# Date 1331532893 0
# Node ID 7ce064380002a1d7ca7f87cd83eaef7fb5b8a725
# Parent  40cc1c0f20a8adcb4fac4a498db08e8f8c5a6c6b
Implement group-delay in CARFAC_Transfer_Functions, and adjust design parameters to make it come out good (with increased min damping in low-f channels)

diff -r 40cc1c0f20a8 -r 7ce064380002 matlab/bmm/carfac/CARFAC_Design.m
--- a/matlab/bmm/carfac/CARFAC_Design.m	Sun Mar 11 22:45:36 2012 +0000
+++ b/matlab/bmm/carfac/CARFAC_Design.m	Mon Mar 12 06:14:53 2012 +0000
@@ -96,7 +96,7 @@
     'v_offset', 0.01, ...  % offset gives a quadratic part
     'v2_corner', 0.2, ...  % corner for essential nonlin
     'v_damp_max', 0.01, ... % damping delta damping from velocity nonlin
-    'min_zeta', 0.10, ...
+    'min_zeta', 0.10, ... % minimum damping factor in mid-freq channels
     'first_pole_theta', 0.85*pi, ...
     'zero_ratio', sqrt(2), ... % how far zero is above pole
     'high_f_damping_compression', 0.5, ... % 0 to 1 to compress zeta
@@ -188,8 +188,12 @@
 %                       and when ff is 1 it goes to zero at theta = pi.
 filter_coeffs.zr_coeffs = zr_coeffs;  % how r relates to zeta
 
-r = (1 - zr_coeffs * filter_params.min_zeta);
-filter_coeffs.r1_coeffs = r;
+min_zeta = filter_params.min_zeta;
+% increase the min damping where channels are spaced out more:
+min_zeta = min_zeta + 0.25*(ERB_Hz(pole_freqs) ./ pole_freqs - min_zeta);
+r1 = (1 - zr_coeffs .* min_zeta);  % "1" for the min-damping condition
+
+filter_coeffs.r1_coeffs = r1;
 
 % undamped coupled-form coefficients:
 filter_coeffs.a0_coeffs = a0;
@@ -200,7 +204,7 @@
 filter_coeffs.h_coeffs = h;
 
 % for unity gain at min damping, radius r; only used in CARFAC_Init:
-extra_damping = zeros(size(r));
+extra_damping = zeros(size(r1));
 % this function needs to take filter_coeffs even if we haven't finished
 % constucting it by putting in the g0_coeffs:
 filter_coeffs.g0_coeffs = CARFAC_Stage_g(filter_coeffs, extra_damping);
@@ -407,7 +411,7 @@
 %                       v_offset: 0.0100
 %                      v2_corner: 0.2000
 %                     v_damp_max: 0.0100
-%                       min_zeta: 0.1200
+%                       min_zeta: 0.1000
 %               first_pole_theta: 2.6704
 %                     zero_ratio: 1.4142
 %     high_f_damping_compression: 0.5000
diff -r 40cc1c0f20a8 -r 7ce064380002 matlab/bmm/carfac/CARFAC_Transfer_Functions.m
--- a/matlab/bmm/carfac/CARFAC_Transfer_Functions.m	Sun Mar 11 22:45:36 2012 +0000
+++ b/matlab/bmm/carfac/CARFAC_Transfer_Functions.m	Mon Mar 12 06:14:53 2012 +0000
@@ -18,11 +18,12 @@
 % limitations under the License.
 
 function [complex_transfns_freqs, ...
-  stage_numerators, stage_denominators] = CARFAC_Transfer_Functions( ...
-  CF, freqs, to_channels, from_channels)
+  stage_numerators, stage_denominators, group_delays] = ...
+  CARFAC_Transfer_Functions(CF, freqs, to_channels, from_channels)
 % function [complex_transfns_freqs, ...
-%   stage_numerators, stage_denominators] = CARFAC_Transfer_Functions( ...
-%   CF, freqs, to_channels, from_channels)
+%   stage_numerators, stage_denominators, group_delays] = ...
+%   CARFAC_Transfer_Functions(CF, freqs, to_channels, from_channels)
+%
 % Return transfer functions as polynomials in z (nums & denoms);
 % And evaluate them at freqs if it's given, to selected output,
 %   optionally from selected starting points (from 0, input, by default).
@@ -49,7 +50,7 @@
   
   % Now multiply gains from input to output places; use logs?
   log_gains = log(gains);
-  cum_log_gains = cumsum(log_gains);
+  cum_log_gains = cumsum(log_gains);  % accum across cascaded stages  
   
   % And figure out which cascade products we want:
   n_ch = CF.n_ch;
@@ -71,6 +72,13 @@
   from_cum(not_input, :) = cum_log_gains(from_channels(not_input), :);
   log_transfns = cum_log_gains(to_channels, :) - from_cum;
   complex_transfns_freqs = exp(log_transfns);
+  
+  if nargout >= 4
+    phases = imag(log_gains);  % no wrapping problem on single stages
+    cum_phases = cumsum(phases);  % so no wrapping here either
+    group_delays = -diff(cum_phases')';  % diff across frequencies
+    group_delays = group_delays ./ (2*pi*repmat(diff(freqs), n_ch, 1));
+  end
 else
   % If no freqs are provided, do nothing but return the stage info above:
   complex_transfns_freqs = [];