dicklyon@615: % Copyright 2013, Google, Inc.
dicklyon@615: % Author: Richard F. Lyon
dicklyon@615: %
dicklyon@615: % This Matlab file is part of an implementation of Lyon's cochlear model:
dicklyon@615: % "Cascade of Asymmetric Resonators with Fast-Acting Compression"
dicklyon@615: % to supplement Lyon's upcoming book "Human and Machine Hearing"
dicklyon@615: %
dicklyon@615: % Licensed under the Apache License, Version 2.0 (the "License");
dicklyon@615: % you may not use this file except in compliance with the License.
dicklyon@615: % You may obtain a copy of the License at
dicklyon@615: %
dicklyon@615: %     http://www.apache.org/licenses/LICENSE-2.0
dicklyon@615: %
dicklyon@615: % Unless required by applicable law or agreed to in writing, software
dicklyon@615: % distributed under the License is distributed on an "AS IS" BASIS,
dicklyon@615: % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
dicklyon@615: % See the License for the specific language governing permissions and
dicklyon@615: % limitations under the License.
dicklyon@615: 
dicklyon@615: function [frame_rate, num_frames] = SAI_RunLayered(CF, input_waves)
ronw@675: % function [CF, SAI_movie] = CARFAC_RunLayered(CF, input_waves)
dicklyon@615: % This function runs the CARFAC and generates an SAI movie, dumped as PNG
dicklyon@615: % files for now.
ronw@676: %
ronw@676: % Computes a "layered" SAI composed of images computed at several
ronw@676: % time scales.
ronw@676: %
dicklyon@615: % Layer 1 is not decimated from the 22050 rate; subsequent layers have
ronw@675: % smoothing and 2X decimation each.  All layers get composited together
dicklyon@615: % into movie frames.
dicklyon@615: 
dicklyon@615: n_ch = CF.n_ch;
dicklyon@615: [n_samp, n_ears] = size(input_waves);
dicklyon@615: if n_ears ~= CF.n_ears
dicklyon@615:   error('bad number of input_waves channels passed to CARFAC_Run')
dicklyon@615: end
dicklyon@615: fs = CF.fs;
dicklyon@615: 
dicklyon@665: seglen = round(fs / 30);  % Pick about 30 fps
dicklyon@665: frame_rate = fs / seglen;
ronw@676: n_segs = ceil(n_samp / seglen);
ronw@676: 
dicklyon@665: 
dicklyon@615: % Design the composite log-lag SAI using these parameters and defaults.
dicklyon@665: n_layers = 15
dicklyon@665: width_per_layer = 36;
dicklyon@665: [layer_array, total_width, lags] = ...
dicklyon@665:   SAI_DesignLayers(n_layers, width_per_layer, seglen);
dicklyon@665: 
dicklyon@665: % Find where in the lag curve corresponds to the piano black keys:
dicklyon@665: pitches = fs ./ lags;
dicklyon@665: key_indices = [];
dicklyon@665: df = log(2)/width_per_layer;
dicklyon@665: for f = [BlackKeyFrequencies, 8, 4, 2, 1-df, 1, 1+df, 0.5, 0.25, 0.125, ...
dicklyon@665:     -2000, -1000, -500, -250, -125];  % Augment with beat.
dicklyon@665:   [dist, index] = min((f - pitches).^2);
dicklyon@665:   key_indices = [key_indices, index];
dicklyon@665: end
dicklyon@665: piano = zeros(1, total_width);
dicklyon@665: piano(key_indices) = 1;
dicklyon@665: piano = [piano; piano; piano];
dicklyon@665: 
dicklyon@615: 
dicklyon@615: % Make the composite SAI image array.
dicklyon@615: composite_frame = zeros(n_ch, total_width);
dicklyon@615: 
dicklyon@615: % Make the history buffers in the layers_array:
dicklyon@615: for layer = 1:n_layers
dicklyon@615:   layer_array(layer).nap_buffer = zeros(layer_array(layer).buffer_width, n_ch);
dicklyon@615:   layer_array(layer).nap_fraction = 0;  % leftover fraction to shift in.
dicklyon@619:   % The SAI frame is transposed to be image-like.
dicklyon@619:   layer_array(layer).frame = zeros(n_ch, layer_array(layer).frame_width);
dicklyon@615: end
dicklyon@615: 
dicklyon@665: n_marginal_rows = 100;
dicklyon@619: marginals = [];
dicklyon@665: average_composite = 0;
dicklyon@619: 
dicklyon@665: future_lags = layer_array(1).future_lags;
dicklyon@665: % marginals_frame = zeros(total_width - future_lags + 2*n_ch, total_width);
dicklyon@665: marginals_frame = zeros(n_ch, total_width);
dicklyon@665: 
dicklyon@615: for seg_num = 1:n_segs
ronw@676:   % seg_range is the range of input sample indices for this segment
dicklyon@615:   if seg_num == n_segs
dicklyon@615:     % The last segment may be short of seglen, but do it anyway:
ronw@676:     seg_range = (seglen*(seg_num - 1) + 1):n_samp;
dicklyon@615:   else
ronw@676:     seg_range = seglen*(seg_num - 1) + (1:seglen);
dicklyon@615:   end
ronw@676:   [seg_naps, CF] = CARFAC_Run_Segment(CF, input_waves(seg_range, :));
dicklyon@615:   
dicklyon@615:   seg_naps = max(0, seg_naps);  % Rectify
dicklyon@615:   
dicklyon@615:   if seg_num == n_segs  % pad out the last result
dicklyon@615:     seg_naps = [seg_naps; zeros(seglen - size(seg_naps,1), size(seg_naps, 2))];
dicklyon@615:   end
dicklyon@615:  
dicklyon@615:   % Shift new data into some or all of the layer buffers:
dicklyon@615:   layer_array = SAI_UpdateBuffers(layer_array, seg_naps, seg_num);
dicklyon@615: 
dicklyon@665:   for layer = n_layers:-1:1  % Stabilize and blend from coarse to fine
dicklyon@615:     update_interval = layer_array(layer).update_interval;
dicklyon@615:     if 0 == mod(seg_num, update_interval)
dicklyon@619:       layer_array(layer) = SAI_StabilizeLayer(layer_array(layer));
dicklyon@619:       composite_frame = SAI_BlendFrameIntoComposite( ...
dicklyon@615:         layer_array(layer), composite_frame);
dicklyon@615:     end
dicklyon@615:   end
dicklyon@665:   
dicklyon@665:   average_composite = average_composite + ...
dicklyon@665:     0.01 * (composite_frame - average_composite);
dicklyon@619:  
dicklyon@619:   if isempty(marginals)
dicklyon@665:     marginals = zeros(n_marginal_rows, total_width);
dicklyon@619:   end
dicklyon@619:   for row = n_marginal_rows:-1:11
dicklyon@619:     % smooth from row above (lower number)
dicklyon@619:     marginals(row, :) = marginals(row, :) + ...
dicklyon@665:       2^((10 - row)/8) * (1.01*marginals(row - 1, :) - marginals(row, :));
dicklyon@619:   end
dicklyon@615:   lag_marginal = mean(composite_frame, 1);  % means max out near 1 or 2
dicklyon@665:   lag_marginal = lag_marginal - 0.75*smooth1d(lag_marginal, 30)';
dicklyon@665:   
dicklyon@665:   freq_marginal = mean(layer_array(1).nap_buffer);
dicklyon@665:   % emphasize local peaks:
dicklyon@665:   freq_marginal = freq_marginal - 0.5*smooth1d(freq_marginal, 5)';
dicklyon@665:   
dicklyon@665:   
dicklyon@665: %   marginals_frame = [marginals_frame(:, 2:end), ...
dicklyon@665: %     [lag_marginal(1:(end - future_lags)), freq_marginal(ceil((1:(2*end))/2))]'];
dicklyon@665:   marginals_frame = [marginals_frame(:, 2:end), freq_marginal(1:end)'];
dicklyon@665:   
dicklyon@619:   for row = 10:-1:1
dicklyon@665:     marginals(row, :) = lag_marginal - (10 - row) / 40;
dicklyon@615:   end
dicklyon@615:     
dicklyon@615:   if 0 == mod(seg_num, update_interval) || seg_num == 1
dicklyon@615:     coc_gram = layer_array(end).nap_buffer';
dicklyon@615:     [n_ch, n_width] = size(composite_frame);
dicklyon@615:     coc_gram = [coc_gram, zeros(n_ch, n_width - size(coc_gram, 2))];
dicklyon@665:     coc_gram = coc_gram(:, (end-total_width+1):end);
dicklyon@615:   end
dicklyon@615:   
dicklyon@665:   display_frame = [ ...  % coc_gram; ...
dicklyon@665:     4 * marginals_frame; ...
dicklyon@665:     composite_frame(ceil((1:(2*end))/2), :); ...
dicklyon@665:     piano; ...
dicklyon@665:     10*max(0,marginals)];
dicklyon@615:   
dicklyon@615:   cmap = jet;
dicklyon@615:   cmap = 1 - gray;  % jet
dicklyon@615:   figure(10)
dicklyon@615:   image(32*display_frame);
dicklyon@615:   colormap(cmap);
dicklyon@615: 
dicklyon@615:   drawnow
dicklyon@615:   imwrite(32*display_frame, cmap, sprintf('frames/frame%05d.png', seg_num));
dicklyon@615: end
dicklyon@615: 
dicklyon@615: num_frames = seg_num;
dicklyon@615: 
dicklyon@615: return
dicklyon@615: 
dicklyon@615: 
dicklyon@665: function frequencies = BlackKeyFrequencies
dicklyon@665: black_indices = [];
dicklyon@665: for index = 0:87
dicklyon@665:   if any(mod(index, 12) == [1 4 6 9 11])
dicklyon@665:     black_indices = [black_indices, index];
dicklyon@665:   end
dicklyon@665: end
dicklyon@665: frequencies = 27.5 * 2.^(black_indices / 12);
dicklyon@615: 
dicklyon@615: