comparison matlab/bmm/carfac/SAI_DesignLayers.m @ 604:ec3a1c74ec54

Add files for making log-lag SAI from CARFAC's NAP output.  The file
SAI_RunLayered.m dumps frames to PNG files.  The hacking script calls
ffmpeg to assemble them with the soundtrack into a movie.

+% Copyright 2013, Google, Inc.
+% Author: Richard F. Lyon
+%
+% This Matlab file is part of an implementation of Lyon's cochlear model:
+% "Cascade of Asymmetric Resonators with Fast-Acting Compression"
+% to supplement Lyon's upcoming book "Human and Machine Hearing"
+%
+% Licensed under the Apache License, Version 2.0 (the "License");
+% you may not use this file except in compliance with the License.
+% You may obtain a copy of the License at
+%
+%     http://www.apache.org/licenses/LICENSE-2.0
+%
+% Unless required by applicable law or agreed to in writing, software
+% distributed under the License is distributed on an "AS IS" BASIS,
+% 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 [layer_array, total_width] = SAI_DesignLayers( ...
+  n_layers, width_per_layer)
+% function [layer_array, total_width] = SAI_DesignLayers( ...
+%   n_layers, width_per_layer)
+%
+% The layer_array is a struct array containing an entry for each layer
+% in a layer of power-of-2 decimated pieces of SAI that get composited
+% into a log-lag SAI.
+% Each struct has the following fields:
+%  .width - number of pixels occupied in the final composite SAI,
+%     not counting the overlap into pixels counted for other layers.
+%  .target_indices - column indices in the final composite SAI,
+%     counting the overlap region(s).
+%  .lag_curve - for each point in the final composite SAI, the float index
+%     in the layer's buffer to interp from.
+%  .alpha - the blending coefficent, mostly 1, tapering toward 0 in the overlap
+%     region(s).
+% Layer 1 has no overlap to it right, and layer n_layers has none to its
+% left, but sizes of the target_indices, lag_curve, and alpha vectors are
+% otherwise width + left_overlap + right_overlap.  The total width of the
+% final composite SAI is the sum of the widths.
+% Other fields could be added to hold state, such as history buffers for
+% each layer, or those could go in state struct array...
+
+
+% Elevate these to a param struct?
+if nargin < 1
+  n_layers = 11
+end
+if nargin < 2
+  width_per_layer = 32;  % resolution "half life" in space
+end
+future_lags = 3 * width_per_layer;
+width_first_layer = future_lags + 2 * width_per_layer;
+width_extra_last_layer = 2 * width_per_layer;
+left_overlap = 15;
+right_overlap = 15;
+first_window_width = 400;  % or maybe use seglen?  or 0.020 * fs?
+min_window_width = 2*width_per_layer;  % or somewhere on that order
+window_exponent = 1.4;
+alpha_max = 0.5;
+
+% Start with NAP_samples_per_SAI_sample, declining to 1 from here:
+max_samples_per = 2^(n_layers - 1);
+% Construct the overall lag-warping function:
+NAP_samples_per_SAI_sample = [ ...
+  max_samples_per * ones(1, width_extra_last_layer), ...
+  max_samples_per * ...
+    2 .^ (-(1:(width_per_layer * (n_layers - 1))) / width_per_layer), ...
+  ones(1, width_first_layer)];
+
+% Each layer needs a lag_warp for a portion of that, divided by
+% 2^(layer-1), where the portion includes some overlap into its neighbors
+% with higher layer numbers on left, lower on right.
+
+% Layer 1, rightmost, representing recent, current and near-future (negative
+% lag) relative to trigger time, has 1 NAP sample per SAI sample.  Other
+% layers map more than one NAP sample into 1 SAI sample.  Layer 2 is
+% computed as 2X decimated, 2 NAP samples per SAI sample, but then gets 
+% interpolated to between 1 and 2 (and outside that range in the overlap
+% regions) to connect up smoothly.  Each layer is another 2X decimated.
+% The last layer limits out at 1 (representing 2^(n_layers) SAI samples)
+% at the width_extra_last_layer SAI samples that extend to the far past.
+
+layer_array = [];  % to hold a struct array
+for layer = 1:n_layers
+  layer_array(layer).width = width_per_layer;
+  layer_array(layer).left_overlap = left_overlap;
+  layer_array(layer).right_overlap = right_overlap;
+  layer_array(layer).future_lags = 0;
+  % Layer decimation factors:  1 1 1 1 2 2 2 4 4 4 8 ...
+  layer_array(layer).update_interval = max(1, 2 ^ floor((layer - 2) / 3));
+end
+% Patch up the exceptions.
+layer_array(1).width = width_first_layer;
+layer_array(end).width = layer_array(end).width + width_extra_last_layer;
+layer_array(1).right_overlap = 0;
+layer_array(end).left_overlap = 0;
+layer_array(1).future_lags = future_lags;
+
+% For each layer, working backwards, from left, find the locations they
+% they render into in the final SAI.
+offset = 0;
+for layer = n_layers:-1:1
+  width = layer_array(layer).width;
+  left = layer_array(layer).left_overlap;
+  right = layer_array(layer).right_overlap;
+  
+  % Size of the vectors needed.
+  n_final_lags = left + width + right;
+  layer_array(layer).n_final_lags = n_final_lags;
+  
+  % Integer indices into the final composite SAI for this layer.
+  target_indices = ((1 - left):(width + right)) + offset;
+  layer_array(layer).target_indices = target_indices;
+    
+  % Make a blending coefficient alpha, ramped in the overlap zone.
+  alpha = ones(1, n_final_lags);
+  alpha(1:left) = alpha(1:left) .* (1:left)/(left + 1);
+  alpha(end + 1 - (1:right)) = ...
+    alpha(end + 1 - (1:right)) .* (1:right)/(right + 1);
+  layer_array(layer).alpha = alpha * alpha_max;
+  
+  offset = offset + width;  % total width from left through this layer.
+end
+total_width = offset;  % Return size of SAI this will make.
+
+% for each layer, fill in its lag-resampling function for interp1:
+for layer = 1:n_layers
+  width = layer_array(layer).width;
+  left = layer_array(layer).left_overlap;
+  right = layer_array(layer).right_overlap;
+  
+  % Still need to adjust this to make lags match at edges:
+  target_indices = layer_array(layer).target_indices;
+  samples_per = NAP_samples_per_SAI_sample(target_indices);
+  % Accumulate lag backwards from the zero-lag point, convert to units of
+  % samples in the current layer.
+  lag_curve = (cumsum(samples_per(end:-1:1))) / 2^(layer-1);
+  lag_curve = lag_curve(end:-1:1);  % Turn it back to corrent order.
+  % Now adjust it to match the zero-lag point or a lag-point from
+  % previous layer, and reverse it back into place.
+  if layer == 1
+    lag_adjust = lag_curve(end) - 0;
+  else
+    % Align right edge to previous layer's left edge, adjusting for 2X
+    % scaling factor difference.
+    lag_adjust = lag_curve(end - right) - last_left_lag / 2;
+  end
+  lag_curve = lag_curve - lag_adjust;
+  % lag_curve is now offsets from right end of layer's frame.
+  layer_array(layer).lag_curve = lag_curve;
+  % Specify number of point to generate in pre-warp frame.
+  layer_array(layer).frame_width = ceil(1 + lag_curve(1));
+  if layer < n_layers  % to avoid the left = 0 unused end case.
+    % A point to align next layer to.
+    last_left_lag = lag_curve(left) - layer_array(layer).future_lags;  
+  end
+  
+  % Specify a good window width (in history buffer, for picking triggers) 
+  % in samples for this layer, exponentially approaching minimum.
+  layer_array(layer).window_width = round(min_window_width + ...
+    first_window_width / window_exponent^(layer - 1));
+  
+  % Say about how long the history buffer needs to be to shift any trigger
+  % location in the range of the window to a fixed location.  Assume
+  % using two window placements overlapped 50%.
+  n_triggers = 2;
+  layer_array(layer).buffer_width = layer_array(layer).frame_width + ...
+    ceil((1 + (n_triggers - 1)/2) * layer_array(layer).window_width);
+end
+
+return
+