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@623: seglen = round(fs / 30);  % Pick about 30 fps
dicklyon@623: frame_rate = fs / seglen;
dicklyon@623: 
dicklyon@604: % Design the composite log-lag SAI using these parameters and defaults.
dicklyon@623: n_layers = 15
dicklyon@623: width_per_layer = 36;
dicklyon@623: [layer_array, total_width, lags] = ...
dicklyon@623:   SAI_DesignLayers(n_layers, width_per_layer, seglen);
dicklyon@623: 
dicklyon@623: % Find where in the lag curve corresponds to the piano black keys:
dicklyon@623: pitches = fs ./ lags;
dicklyon@623: key_indices = [];
dicklyon@623: df = log(2)/width_per_layer;
dicklyon@623: for f = [BlackKeyFrequencies, 8, 4, 2, 1-df, 1, 1+df, 0.5, 0.25, 0.125, ...
dicklyon@623:     -2000, -1000, -500, -250, -125];  % Augment with beat.
dicklyon@623:   [dist, index] = min((f - pitches).^2);
dicklyon@623:   key_indices = [key_indices, index];
dicklyon@623: end
dicklyon@623: piano = zeros(1, total_width);
dicklyon@623: piano(key_indices) = 1;
dicklyon@623: piano = [piano; piano; piano];
dicklyon@623: 
dicklyon@604: 
dicklyon@604: % Make the composite SAI image array.
dicklyon@604: composite_frame = zeros(n_ch, total_width);
dicklyon@604: 
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@608:   % The SAI frame is transposed to be image-like.
dicklyon@608:   layer_array(layer).frame = zeros(n_ch, layer_array(layer).frame_width);
dicklyon@604: end
dicklyon@604: 
dicklyon@623: n_marginal_rows = 100;
dicklyon@608: marginals = [];
dicklyon@623: average_composite = 0;
dicklyon@608: 
dicklyon@623: future_lags = layer_array(1).future_lags;
dicklyon@623: % marginals_frame = zeros(total_width - future_lags + 2*n_ch, total_width);
dicklyon@623: marginals_frame = zeros(n_ch, total_width);
dicklyon@623: 
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@623:   for layer = n_layers:-1:1  % Stabilize and blend from coarse to fine
dicklyon@604:     update_interval = layer_array(layer).update_interval;
dicklyon@604:     if 0 == mod(seg_num, update_interval)
dicklyon@608:       layer_array(layer) = SAI_StabilizeLayer(layer_array(layer));
dicklyon@608:       composite_frame = SAI_BlendFrameIntoComposite( ...
dicklyon@604:         layer_array(layer), composite_frame);
dicklyon@604:     end
dicklyon@604:   end
dicklyon@623:   
dicklyon@623:   average_composite = average_composite + ...
dicklyon@623:     0.01 * (composite_frame - average_composite);
dicklyon@608:  
dicklyon@608:   if isempty(marginals)
dicklyon@623:     marginals = zeros(n_marginal_rows, total_width);
dicklyon@608:   end
dicklyon@608:   for row = n_marginal_rows:-1:11
dicklyon@608:     % smooth from row above (lower number)
dicklyon@608:     marginals(row, :) = marginals(row, :) + ...
dicklyon@623:       2^((10 - row)/8) * (1.01*marginals(row - 1, :) - marginals(row, :));
dicklyon@608:   end
dicklyon@604:   lag_marginal = mean(composite_frame, 1);  % means max out near 1 or 2
dicklyon@623:   lag_marginal = lag_marginal - 0.75*smooth1d(lag_marginal, 30)';
dicklyon@623:   
dicklyon@623:   freq_marginal = mean(layer_array(1).nap_buffer);
dicklyon@623:   % emphasize local peaks:
dicklyon@623:   freq_marginal = freq_marginal - 0.5*smooth1d(freq_marginal, 5)';
dicklyon@623:   
dicklyon@623:   
dicklyon@623: %   marginals_frame = [marginals_frame(:, 2:end), ...
dicklyon@623: %     [lag_marginal(1:(end - future_lags)), freq_marginal(ceil((1:(2*end))/2))]'];
dicklyon@623:   marginals_frame = [marginals_frame(:, 2:end), freq_marginal(1:end)'];
dicklyon@623:   
dicklyon@608:   for row = 10:-1:1
dicklyon@623:     marginals(row, :) = lag_marginal - (10 - row) / 40;
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@623:     coc_gram = coc_gram(:, (end-total_width+1):end);
dicklyon@604:   end
dicklyon@604:   
dicklyon@623:   display_frame = [ ...  % coc_gram; ...
dicklyon@623:     4 * marginals_frame; ...
dicklyon@623:     composite_frame(ceil((1:(2*end))/2), :); ...
dicklyon@623:     piano; ...
dicklyon@623:     10*max(0,marginals)];
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@623: function frequencies = BlackKeyFrequencies
dicklyon@623: black_indices = [];
dicklyon@623: for index = 0:87
dicklyon@623:   if any(mod(index, 12) == [1 4 6 9 11])
dicklyon@623:     black_indices = [black_indices, index];
dicklyon@623:   end
dicklyon@623: end
dicklyon@623: frequencies = 27.5 * 2.^(black_indices / 12);
dicklyon@604: 
dicklyon@604: