annotate trunk/matlab/bmm/carfac/SAI_RunLayered.m @ 619:2e456754fe20

Better functionization of SAI, and new way to make picture with lag marginal and smoothed history of lag marginal.
author dicklyon@google.com
date Mon, 13 May 2013 21:15:56 +0000
parents 2b2ef398b557
children d0ff15c36828
rev   line source
dicklyon@615 1 % Copyright 2013, Google, Inc.
dicklyon@615 2 % Author: Richard F. Lyon
dicklyon@615 3 %
dicklyon@615 4 % This Matlab file is part of an implementation of Lyon's cochlear model:
dicklyon@615 5 % "Cascade of Asymmetric Resonators with Fast-Acting Compression"
dicklyon@615 6 % to supplement Lyon's upcoming book "Human and Machine Hearing"
dicklyon@615 7 %
dicklyon@615 8 % Licensed under the Apache License, Version 2.0 (the "License");
dicklyon@615 9 % you may not use this file except in compliance with the License.
dicklyon@615 10 % You may obtain a copy of the License at
dicklyon@615 11 %
dicklyon@615 12 % http://www.apache.org/licenses/LICENSE-2.0
dicklyon@615 13 %
dicklyon@615 14 % Unless required by applicable law or agreed to in writing, software
dicklyon@615 15 % distributed under the License is distributed on an "AS IS" BASIS,
dicklyon@615 16 % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
dicklyon@615 17 % See the License for the specific language governing permissions and
dicklyon@615 18 % limitations under the License.
dicklyon@615 19
dicklyon@615 20 function [frame_rate, num_frames] = SAI_RunLayered(CF, input_waves)
dicklyon@615 21 % function [CF, SAI_movie] = CARFAC_Run_Layered_SAI(CF, input_waves)
dicklyon@615 22 % This function runs the CARFAC and generates an SAI movie, dumped as PNG
dicklyon@615 23 % files for now.
dicklyon@615 24
dicklyon@615 25 % Layer 1 is not decimated from the 22050 rate; subsequent layers have
dicklyon@615 26 % smoothing and 2X decimation each. All layers get composited togehter
dicklyon@615 27 % into movie frames.
dicklyon@615 28
dicklyon@615 29 n_ch = CF.n_ch;
dicklyon@615 30 [n_samp, n_ears] = size(input_waves);
dicklyon@615 31 if n_ears ~= CF.n_ears
dicklyon@615 32 error('bad number of input_waves channels passed to CARFAC_Run')
dicklyon@615 33 end
dicklyon@615 34 fs = CF.fs;
dicklyon@615 35
dicklyon@615 36 % Design the composite log-lag SAI using these parameters and defaults.
dicklyon@615 37 n_layers = 10;
dicklyon@615 38 width_per_layer = 40;
dicklyon@615 39 [layer_array, total_width] = SAI_DesignLayers(n_layers, width_per_layer);
dicklyon@615 40
dicklyon@615 41 % Make the composite SAI image array.
dicklyon@615 42 composite_frame = zeros(n_ch, total_width);
dicklyon@615 43
dicklyon@619 44 seglen = round(fs / 30); % Pick about 60 fps
dicklyon@615 45 frame_rate = fs / seglen;
dicklyon@615 46 n_segs = ceil(n_samp / seglen);
dicklyon@615 47
dicklyon@615 48 % Make the history buffers in the layers_array:
dicklyon@615 49 for layer = 1:n_layers
dicklyon@615 50 layer_array(layer).nap_buffer = zeros(layer_array(layer).buffer_width, n_ch);
dicklyon@615 51 layer_array(layer).nap_fraction = 0; % leftover fraction to shift in.
dicklyon@619 52 % The SAI frame is transposed to be image-like.
dicklyon@619 53 layer_array(layer).frame = zeros(n_ch, layer_array(layer).frame_width);
dicklyon@615 54 end
dicklyon@615 55
dicklyon@619 56 n_marginal_rows = 34;
dicklyon@619 57 marginals = [];
dicklyon@619 58
dicklyon@619 59 average_frame = 0;
dicklyon@615 60 for seg_num = 1:n_segs
dicklyon@615 61 % k_range is the range of input sample indices for this segment
dicklyon@615 62 if seg_num == n_segs
dicklyon@615 63 % The last segment may be short of seglen, but do it anyway:
dicklyon@615 64 k_range = (seglen*(seg_num - 1) + 1):n_samp;
dicklyon@615 65 else
dicklyon@615 66 k_range = seglen*(seg_num - 1) + (1:seglen);
dicklyon@615 67 end
dicklyon@615 68 % Process a segment to get a slice of decim_naps, and plot AGC state:
dicklyon@615 69 [seg_naps, CF] = CARFAC_Run_Segment(CF, input_waves(k_range, :));
dicklyon@615 70
dicklyon@615 71 seg_naps = max(0, seg_naps); % Rectify
dicklyon@615 72
dicklyon@615 73 if seg_num == n_segs % pad out the last result
dicklyon@615 74 seg_naps = [seg_naps; zeros(seglen - size(seg_naps,1), size(seg_naps, 2))];
dicklyon@615 75 end
dicklyon@615 76
dicklyon@615 77 % Shift new data into some or all of the layer buffers:
dicklyon@615 78 layer_array = SAI_UpdateBuffers(layer_array, seg_naps, seg_num);
dicklyon@615 79
dicklyon@615 80
dicklyon@615 81 for layer = n_layers:-1:1 % blend from coarse to fine
dicklyon@615 82 update_interval = layer_array(layer).update_interval;
dicklyon@615 83 if 0 == mod(seg_num, update_interval)
dicklyon@619 84 layer_array(layer) = SAI_StabilizeLayer(layer_array(layer));
dicklyon@619 85 new_frame = layer_array(layer).frame;
dicklyon@619 86 composite_frame = SAI_BlendFrameIntoComposite( ...
dicklyon@615 87 layer_array(layer), composite_frame);
dicklyon@615 88 end
dicklyon@615 89 end
dicklyon@619 90
dicklyon@619 91 if isempty(marginals)
dicklyon@619 92 composite_width = size(composite_frame, 2);
dicklyon@619 93 marginals = zeros(n_marginal_rows, composite_width);
dicklyon@619 94 end
dicklyon@619 95 for row = n_marginal_rows:-1:11
dicklyon@619 96 % smooth from row above (lower number)
dicklyon@619 97 marginals(row, :) = marginals(row, :) + ...
dicklyon@619 98 2^((10 - row)/2) * (1.06*marginals(row - 1, :) - marginals(row, :));
dicklyon@619 99 end
dicklyon@615 100 lag_marginal = mean(composite_frame, 1); % means max out near 1 or 2
dicklyon@619 101 for row = 10:-1:1
dicklyon@619 102 marginals(row, :) = (lag_marginal - smooth1d(lag_marginal, 30)') - ...
dicklyon@619 103 (10 - row) / 40;
dicklyon@615 104 end
dicklyon@615 105
dicklyon@615 106 if 0 == mod(seg_num, update_interval) || seg_num == 1
dicklyon@615 107 coc_gram = layer_array(end).nap_buffer';
dicklyon@615 108 [n_ch, n_width] = size(composite_frame);
dicklyon@615 109 coc_gram = [coc_gram, zeros(n_ch, n_width - size(coc_gram, 2))];
dicklyon@615 110 end
dicklyon@615 111
dicklyon@619 112 display_frame = [coc_gram; ...
dicklyon@619 113 composite_frame(floor(1:0.5:end), :); 20*max(0,marginals)];
dicklyon@615 114
dicklyon@615 115 cmap = jet;
dicklyon@615 116 cmap = 1 - gray; % jet
dicklyon@615 117 figure(10)
dicklyon@615 118 image(32*display_frame);
dicklyon@615 119 colormap(cmap);
dicklyon@615 120
dicklyon@615 121 drawnow
dicklyon@615 122 imwrite(32*display_frame, cmap, sprintf('frames/frame%05d.png', seg_num));
dicklyon@615 123 end
dicklyon@615 124
dicklyon@615 125 num_frames = seg_num;
dicklyon@615 126
dicklyon@615 127 return
dicklyon@615 128
dicklyon@615 129
dicklyon@615 130
dicklyon@615 131