annotate trunk/matlab/bmm/carfac/SAI_StabilizeLayer.m @ 705:be2b68ced23d

Make the C++ SAI match the Matlab implementation.

 - Fix syntax error which caused a linear window to be used instead of the correct sine window.
 - Correct off-by-one indexing errors to compensate for 1-based indexing in Matlab.
 - Enable test to verify that the C++ output matches Matlab with high precision.

function layer_struct = SAI_StabilizeLayer(layer_struct)
% Pick trigger points in buffer, shift rows to offset_from_end,
% and blend into frame

frame = layer_struct.frame;

nap_buffer = real(layer_struct.nap_buffer);
n_buffer_times = size(nap_buffer, 1);
[n_ch, width] = size(frame);

% Make the window to use for all the channels at this layer.
window_width = layer_struct.window_width;
n_window_pos = layer_struct.n_window_pos;
% Windows are always (approx) 50% overlapped:
window_hop = window_width / 2;

window = sin((1:window_width)' * pi / window_width);
window_start = (n_buffer_times - window_width) - ...
    floor((n_window_pos - 1) * window_hop);
window_range = (1:window_width) + window_start - layer_struct.future_lags;
% This should not go negative!
offset_range = (1:width) + window_start - width;
% CHECK
if any(offset_range < 0)
  error;
end

% smooth across channels; more in later layers
smoothed_buffer = smooth1d(nap_buffer', layer_struct.channel_smoothing_scale)';

% For each buffer column (channel), pick a trigger and align into SAI_frame
for ch = 1:n_ch
  smooth_wave = smoothed_buffer(:, ch);  % for the trigger
  
  % Do several window positions and triggers
  for w = 1:n_window_pos
    % move the window to later and go again
    current_window_offset = floor((w - 1) * window_hop);
    [peak_val, trigger_time] = ...
        max(smooth_wave(window_range + current_window_offset) .* window);
    nap_wave = nap_buffer(:, ch);  % for the waveform
    if peak_val <= 0  % just use window center instead
      [peak_val, trigger_time] = max(window);
    end
    trigger_time = trigger_time + current_window_offset;
    alpha = (0.025 + peak_val) / (0.5 + peak_val);  % alpha 0.05 to near 1.0
    frame(ch, :) = alpha * nap_wave(trigger_time + offset_range)' + ...
      (1 - alpha) * frame(ch, :);
  end
end

layer_struct.frame = frame;