tom@516: % Copyright 2012, Google, Inc.
tom@516: % Author: Richard F. Lyon
tom@516: %
tom@516: % This Matlab file is part of an implementation of Lyon's cochlear model:
tom@516: % "Cascade of Asymmetric Resonators with Fast-Acting Compression"
tom@516: % to supplement Lyon's upcoming book "Human and Machine Hearing"
tom@516: %
tom@516: % Licensed under the Apache License, Version 2.0 (the "License");
tom@516: % you may not use this file except in compliance with the License.
tom@516: % You may obtain a copy of the License at
tom@516: %
tom@516: %     http://www.apache.org/licenses/LICENSE-2.0
tom@516: %
tom@516: % Unless required by applicable law or agreed to in writing, software
tom@516: % distributed under the License is distributed on an "AS IS" BASIS,
tom@516: % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
tom@516: % See the License for the specific language governing permissions and
tom@516: % limitations under the License.
tom@516: 
tom@516: function [naps, CF, decim_naps] = CARFAC_Run ...
tom@516:   (CF, input_waves, AGC_plot_fig_num)
tom@516: % function [naps, CF, CF.cum_k, decim_naps] = CARFAC_Run ...
tom@516: %    (CF, input_waves, CF.cum_k, AGC_plot_fig_num)
tom@516: % This function runs the CARFAC; that is, filters a 1 or more channel
tom@516: % sound input to make one or more neural activity patterns (naps).
tom@516: %
tom@516: % The CF struct holds the filterbank design and state; if you want to
tom@516: % break the input up into segments, you need to use the updated CF
tom@516: % to keep the state between segments.
tom@516: %
tom@516: % input_waves is a column vector if there's just one audio channel;
tom@516: % more generally, it has a row per time sample, a column per audio channel.
tom@516: %
tom@516: % naps has a row per time sample, a column per filterbank channel, and
tom@516: % a layer per audio channel if more than 1.
tom@516: % decim_naps is like naps but time-decimated by the int CF.decimation.
tom@516: %
tom@516: % the input_waves are assumed to be sampled at the same rate as the
tom@516: % CARFAC is designed for; a resampling may be needed before calling this.
tom@516: %
tom@516: % The function works as an outer iteration on time, updating all the
tom@516: % filters and AGC states concurrently, so that the different channels can
tom@516: % interact easily.  The inner loops are over filterbank channels, and
tom@516: % this level should be kept efficient.
tom@516: %
tom@516: % See other functions for designing and characterizing the CARFAC:
tom@516: % CF = CARFAC_Design(fs, CF_filter_params, CF_AGC_params, n_mics)
tom@516: % transfns = CARFAC_Transfer_Functions(CF, to_chans, from_chans)
tom@516: 
tom@516: [n_samp, n_mics] = size(input_waves);
tom@516: n_ch = CF.n_ch;
tom@516: 
tom@516: if nargin < 3
tom@516:   AGC_plot_fig_num = 0;
tom@516: end
tom@516: 
tom@516: if n_mics ~= CF.n_mics
tom@516:   error('bad number of input_waves channels passed to CARFAC_Run')
tom@516: end
tom@516: 
tom@516: % pull coeffs out of struct first, into local vars for convenience
tom@516: decim = CF.AGC_params.decimation;
tom@516: 
tom@516: naps = zeros(n_samp, n_ch, n_mics);
tom@516: if nargout > 2
tom@516:   % make decimated detect output:
tom@516:   decim_naps = zeros(ceil(n_samp/decim), CF.n_ch, CF.n_mics);
tom@516: else
tom@516:   decim_naps = [];
tom@516: end
tom@516: 
tom@516: decim_k = 0;
tom@516: 
tom@516: sum_abs_response = 0;
tom@516: 
tom@516: for k = 1:n_samp
tom@516:   CF.k_mod_decim = mod(CF.k_mod_decim + 1, decim);  % global time phase
tom@516:   % at each time step, possibly handle multiple channels
tom@516:   for mic = 1:n_mics
tom@516:     [filters_out, CF.filter_state(mic)] = CARFAC_FilterStep( ...
tom@516:       input_waves(k, mic), CF.filter_coeffs, CF.filter_state(mic));
tom@516: 
tom@516:     % update IHC state & output on every time step, too
tom@516:     [ihc_out, CF.IHC_state(mic)] = CARFAC_IHCStep( ...
tom@516:       filters_out, CF.IHC_coeffs, CF.IHC_state(mic));
tom@516: 
tom@516: %     sum_abs_response = sum_abs_response + abs(filters_out);
tom@516: 
tom@516:     naps(k, :, mic) = ihc_out;  % output to neural activity pattern
tom@516:   end
tom@516: 
tom@516:   % conditionally update all the AGC stages and channels now:
tom@516:   if CF.k_mod_decim == 0
tom@516:     % just for the plotting option:
tom@516:     decim_k = decim_k + 1;   % index of decimated signal for display
tom@516:     if ~isempty(decim_naps)
tom@516:       for mic = 1:n_mics
tom@516:         % this is HWR out of filters, not IHCs
tom@516:         avg_detect = CF.filter_state(mic).detect_accum / decim;
tom@516:         % This HACK is the IHC version:
tom@516:         avg_detect = CF.IHC_state(mic).ihc_accum / decim;  % for cochleagram
tom@516:         decim_naps(decim_k, :, mic) = avg_detect;  % for cochleagram
tom@516: %         decim_naps(decim_k, :, mic) = sum_abs_response / decim;  % HACK for mechanical out ABS
tom@516: %         sum_abs_response(:) = 0;
tom@516:       end
tom@516:     end
tom@516: 
tom@516:     % get the avg_detects to connect filter_state to AGC_state:
tom@516:     avg_detects = zeros(n_ch, n_mics);
tom@516:     for mic = 1:n_mics
tom@516: %       % mechanical response from filter output through HWR as AGC in:
tom@516: %       avg_detects(:, mic) = CF.filter_state(mic).detect_accum / decim;
tom@516:       CF.filter_state(mic).detect_accum(:) = 0;  % zero the detect accumulator
tom@516:       % New HACK, IHC output relative to rest as input to AGC:
tom@516:       avg_detects(:, mic) = CF.IHC_state(mic).ihc_accum / decim;
tom@516:       CF.IHC_state(mic).ihc_accum(:) = 0;  % zero the detect accumulator
tom@516:     end
tom@516: 
tom@516:     % run the AGC update step:
tom@516:     CF.AGC_state = CARFAC_AGCStep(CF.AGC_coeffs, avg_detects, CF.AGC_state);
tom@516: 
tom@516:     % connect the feedback from AGC_state to filter_state:
tom@516:     for mic = 1:n_mics
tom@516:       new_damping = CF.AGC_state(mic).sum_AGC;
tom@516: %       max_damping = 0.15;  % HACK
tom@516: %       new_damping = min(new_damping, max_damping);
tom@516:       % set the delta needed to get to new_damping:
tom@516:       CF.filter_state(mic).dzB_memory = ...
tom@516:         (new_damping - CF.filter_state(mic).zB_memory) ...
tom@516:           / decim;
tom@516:     end
tom@516: 
tom@516:     if AGC_plot_fig_num
tom@516:       figure(AGC_plot_fig_num); hold off
tom@516:       maxsum = 0;
tom@516:       for mic = 1:n_mics
tom@516:         plot(CF.AGC_state(mic).AGC_memory)
tom@516:         agcsum = sum(CF.AGC_state(mic).AGC_memory, 2);
tom@516:         maxsum(mic) = max(maxsum, max(agcsum));
tom@516:         hold on
tom@516:         plot(agcsum, 'k-')
tom@516:       end
tom@516:       axis([0, CF.n_ch, 0, max(0.001, maxsum)]);
tom@516:       drawnow
tom@516:     end
tom@516:   end
tom@516: end
tom@516: