dicklyon@604: % Copyright 2013, Google, Inc. dicklyon@604: % Author: Richard F. Lyon dicklyon@604: % dicklyon@604: % This Matlab file is part of an implementation of Lyon's cochlear model: dicklyon@604: % "Cascade of Asymmetric Resonators with Fast-Acting Compression" dicklyon@604: % to supplement Lyon's upcoming book "Human and Machine Hearing" dicklyon@604: % dicklyon@604: % Licensed under the Apache License, Version 2.0 (the "License"); dicklyon@604: % you may not use this file except in compliance with the License. dicklyon@604: % You may obtain a copy of the License at dicklyon@604: % dicklyon@604: % http://www.apache.org/licenses/LICENSE-2.0 dicklyon@604: % dicklyon@604: % Unless required by applicable law or agreed to in writing, software dicklyon@604: % distributed under the License is distributed on an "AS IS" BASIS, dicklyon@604: % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. dicklyon@604: % See the License for the specific language governing permissions and dicklyon@604: % limitations under the License. dicklyon@604: dicklyon@604: function [frame_rate, num_frames] = SAI_RunLayered(CF, input_waves) dicklyon@604: % function [CF, SAI_movie] = CARFAC_Run_Layered_SAI(CF, input_waves) dicklyon@604: % This function runs the CARFAC and generates an SAI movie, dumped as PNG dicklyon@604: % files for now. dicklyon@604: dicklyon@604: % Layer 1 is not decimated from the 22050 rate; subsequent layers have dicklyon@604: % smoothing and 2X decimation each. All layers get composited togehter dicklyon@604: % into movie frames. dicklyon@604: dicklyon@604: n_ch = CF.n_ch; dicklyon@604: [n_samp, n_ears] = size(input_waves); dicklyon@604: if n_ears ~= CF.n_ears dicklyon@604: error('bad number of input_waves channels passed to CARFAC_Run') dicklyon@604: end dicklyon@604: fs = CF.fs; dicklyon@604: dicklyon@604: % Design the composite log-lag SAI using these parameters and defaults. dicklyon@604: n_layers = 10; dicklyon@604: width_per_layer = 40; dicklyon@604: [layer_array, total_width] = SAI_DesignLayers(n_layers, width_per_layer); dicklyon@604: dicklyon@604: % Make the composite SAI image array. dicklyon@604: composite_frame = zeros(n_ch, total_width); dicklyon@604: dicklyon@604: seglen = round(fs * 0.020); % Pick about 20 ms segments dicklyon@604: frame_rate = fs / seglen; dicklyon@604: n_segs = ceil(n_samp / seglen); dicklyon@604: dicklyon@604: % Make the history buffers in the layers_array: dicklyon@604: for layer = 1:n_layers dicklyon@604: layer_array(layer).nap_buffer = zeros(layer_array(layer).buffer_width, n_ch); dicklyon@604: layer_array(layer).nap_fraction = 0; % leftover fraction to shift in. dicklyon@604: end dicklyon@604: dicklyon@604: for seg_num = 1:n_segs dicklyon@604: % k_range is the range of input sample indices for this segment dicklyon@604: if seg_num == n_segs dicklyon@604: % The last segment may be short of seglen, but do it anyway: dicklyon@604: k_range = (seglen*(seg_num - 1) + 1):n_samp; dicklyon@604: else dicklyon@604: k_range = seglen*(seg_num - 1) + (1:seglen); dicklyon@604: end dicklyon@604: % Process a segment to get a slice of decim_naps, and plot AGC state: dicklyon@604: [seg_naps, CF] = CARFAC_Run_Segment(CF, input_waves(k_range, :)); dicklyon@604: dicklyon@604: seg_naps = max(0, seg_naps); % Rectify dicklyon@604: dicklyon@604: if seg_num == n_segs % pad out the last result dicklyon@604: seg_naps = [seg_naps; zeros(seglen - size(seg_naps,1), size(seg_naps, 2))]; dicklyon@604: end dicklyon@604: dicklyon@604: % Shift new data into some or all of the layer buffers: dicklyon@604: layer_array = SAI_UpdateBuffers(layer_array, seg_naps, seg_num); dicklyon@604: dicklyon@604: dicklyon@604: for layer = n_layers:-1:1 % blend from coarse to fine dicklyon@604: update_interval = layer_array(layer).update_interval; dicklyon@604: if 0 == mod(seg_num, update_interval) dicklyon@604: nap_buffer = real(layer_array(layer).nap_buffer); dicklyon@604: n_buffer_times = size(nap_buffer, 1); dicklyon@604: width = layer_array(layer).frame_width; % To render linear SAI to. dicklyon@604: new_frame = zeros(n_ch, width); dicklyon@604: dicklyon@604: % Make the window to use for all the channels at this layer. dicklyon@604: layer_factor = 1.5; dicklyon@604: window_size = layer_array(layer).window_width; dicklyon@604: after_samples = layer_array(layer).future_lags; dicklyon@604: dicklyon@604: window_range = (1:window_size) + ... dicklyon@604: (n_buffer_times - window_size) - after_samples; dicklyon@604: window = sin((1:window_size)' * pi / window_size); dicklyon@604: % This should not go negative! dicklyon@604: offset_range = (1:width) + ... dicklyon@604: (n_buffer_times - width - window_size); dicklyon@604: % CHECK dicklyon@604: if any(offset_range < 0) dicklyon@604: error; dicklyon@604: end dicklyon@604: dicklyon@604: % smooth across channels; more in later layers dicklyon@604: smoothed_buffer = smooth1d(nap_buffer', 0.25*(layer - 2))'; dicklyon@604: dicklyon@604: % For each buffer column (channel), pick a trigger and align into SAI_frame dicklyon@604: for ch = 1:n_ch dicklyon@604: smooth_wave = smoothed_buffer(:, ch); % for the trigger dicklyon@604: dicklyon@604: [peak_val, trigger_time] = max(smooth_wave(window_range) .* window); dicklyon@604: nap_wave = nap_buffer(:, ch); % for the waveform dicklyon@604: if peak_val <= 0 % just use window center instead dicklyon@604: [peak_val, trigger_time] = max(window); dicklyon@604: end dicklyon@604: if layer == n_layers % mark the trigger points to display as imaginary. dicklyon@604: layer_array(layer).nap_buffer(trigger_time + window_range(1) - 1, ch) = ... dicklyon@604: layer_array(layer).nap_buffer(trigger_time + window_range(1) - 1, ch) + 1i; dicklyon@604: end dicklyon@604: new_frame(ch, :) = nap_wave(trigger_time + offset_range)'; dicklyon@604: end dicklyon@604: composite_frame = SAI_BlendFrameIntoComposite(new_frame, ... dicklyon@604: layer_array(layer), composite_frame); dicklyon@604: end dicklyon@604: end dicklyon@604: dicklyon@604: dicklyon@604: lag_marginal = mean(composite_frame, 1); % means max out near 1 or 2 dicklyon@604: frame_bottom = zeros(size(composite_frame)); % will end up being 1/3 dicklyon@604: n_bottom_rows = size(frame_bottom, 1); dicklyon@604: for height = 1:n_bottom_rows dicklyon@604: big_ones = lag_marginal > 1*height/n_bottom_rows; dicklyon@604: frame_bottom(n_bottom_rows - height + 1, big_ones) = 2; % 2 for black dicklyon@604: end dicklyon@604: dicklyon@604: if 0 == mod(seg_num, update_interval) || seg_num == 1 dicklyon@604: coc_gram = layer_array(end).nap_buffer'; dicklyon@604: [n_ch, n_width] = size(composite_frame); dicklyon@604: coc_gram = [coc_gram, zeros(n_ch, n_width - size(coc_gram, 2))]; dicklyon@604: trigger_gram = 2 * (imag(coc_gram) ~= 0); dicklyon@604: coc_gram = real(coc_gram); dicklyon@604: end dicklyon@604: dicklyon@604: display_frame = [coc_gram; trigger_gram; ... dicklyon@604: composite_frame(floor(1:0.5:end), :); frame_bottom]; dicklyon@604: dicklyon@604: cmap = jet; dicklyon@604: cmap = 1 - gray; % jet dicklyon@604: figure(10) dicklyon@604: image(32*display_frame); dicklyon@604: colormap(cmap); dicklyon@604: dicklyon@604: drawnow dicklyon@604: imwrite(32*display_frame, cmap, sprintf('frames/frame%05d.png', seg_num)); dicklyon@604: end dicklyon@604: dicklyon@604: num_frames = seg_num; dicklyon@604: dicklyon@604: return dicklyon@604: dicklyon@604: dicklyon@604: dicklyon@604: