comparison trunk/matlab/bmm/carfac/CARFAC_Transfer_Functions.m @ 533:55c46c01e522

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)

 % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 % See the License for the specific language governing permissions and
 % 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).
 %   complex_transfns_freqs has a row of complex gains per to_channel.
 

   z_row = exp((i * 2 * pi / CF.fs) * freqs);  % z = exp(sT)
   gains = Rational_Eval(stage_numerators, stage_denominators, z_row);
   
   % 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;
   if nargin < 3
     to_channels = 1:n_ch;

   from_cum = zeros(length(to_channels), length(z_row));
   not_input = from_channels > 0;
   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 = [];
 end