Mercurial > hg > aimc

changeset 536:2964a3b4a00a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
New design params, including narrower AGC, Greenwood map for more channels, default 71, some renaming, open loop feature, ...
author dicklyon@google.com
date Thu, 22 Mar 2012 22:37:56 +0000
parents dcf18d03d608
children 43c5664a00ec
files trunk/matlab/bmm/carfac/CARFAC_AGC_Step.m trunk/matlab/bmm/carfac/CARFAC_CAR_Step.m trunk/matlab/bmm/carfac/CARFAC_Design.m trunk/matlab/bmm/carfac/CARFAC_Init.m trunk/matlab/bmm/carfac/CARFAC_Run.m trunk/matlab/bmm/carfac/CARFAC_Run_Linear.m trunk/matlab/bmm/carfac/CARFAC_Run_Open_Loop.m trunk/matlab/bmm/carfac/CARFAC_Run_Segment.m trunk/matlab/bmm/carfac/CARFAC_Spatial_Smooth.m trunk/matlab/bmm/carfac/CARFAC_hacking.m trunk/matlab/bmm/carfac/MultiScaleSmooth.m
diffstat 11 files changed, 209 insertions(+), 77 deletions(-) [+]
line wrap: on
line diff
--- a/trunk/matlab/bmm/carfac/CARFAC_AGC_Step.m	Fri Mar 16 04:31:56 2012 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_AGC_Step.m	Thu Mar 22 22:37:56 2012 +0000
@@ -71,9 +71,7 @@
   
   for ear = 1:n_ears
     AGC_in = ins(:,ear);  % the newly decimated input for this ear
-    
-%     AGC_in = max(0, AGC_in);  % don't let neg inputs in
-    
+        
     %  add the latest output (state) of next stage...
     if stage < n_stages
       AGC_in = AGC_in + stage_gain * state(ear).AGC_memory(:, stage+1);
--- a/trunk/matlab/bmm/carfac/CARFAC_CAR_Step.m	Fri Mar 16 04:31:56 2012 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_CAR_Step.m	Thu Mar 22 22:37:56 2012 +0000
@@ -71,7 +71,3 @@
 state.zY_memory = zY;
 state.g_memory = g;
 
-% accum the straight hwr version in case someone wants it:
-hwr_detect = max(0, zY);  % detect with HWR
-state.detect_accum = state.detect_accum + hwr_detect;
-
--- a/trunk/matlab/bmm/carfac/CARFAC_Design.m	Fri Mar 16 04:31:56 2012 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_Design.m	Thu Mar 22 22:37:56 2012 +0000
@@ -74,7 +74,8 @@
   % ERB = 24.7 * (1 + 4.37 * CF_Hz / 1000);
   ERB_Q = 1000/(24.7*4.37);  % 9.2645
   if nargin < 4
-    ERB_break_freq = 1000/4.37;  % 228.833
+%     ERB_break_freq = 1000/4.37;  % 228.833 G&M
+    ERB_break_freq = 165.3;  % Greenwood map's break freq.
   end
 end
 
@@ -84,9 +85,9 @@
     'time_constants', [1, 4, 16, 64]*0.002, ...
     'AGC_stage_gain', 2, ...  % gain from each stage to next slower stage
     'decimation', [8, 2, 2, 2], ...  % how often to update the AGC states
-    'AGC1_scales', [1, 2, 4, 6]*1, ...   % in units of channels
-    'AGC2_scales', [1, 2, 4, 6]*1.5, ... % spread more toward base
-    'detect_scale', 0.15, ...  % the desired damping range
+    'AGC1_scales', [1.0, 1.4,  2.0, 2.8], ...   % in units of channels
+    'AGC2_scales', [1.6, 2.25, 3.2, 4.5], ... % spread more toward base
+    'detect_scale', 0.25, ...  % the desired damping range
     'AGC_mix_coeff', 0.5);
 end
 
@@ -285,7 +286,7 @@
   % when FIR_OK, store the resulting FIR design in coeffs:
   AGC_coeffs.AGC_spatial_iterations(stage) = n_iterations;
   AGC_coeffs.AGC_spatial_FIR(:,stage) = AGC_spatial_FIR;
-  AGC_coeffs.AGC_n_taps(stage) = n_taps;
+  AGC_coeffs.AGC_spatial_n_taps(stage) = n_taps;
   
   % accumulate DC gains from all the stages, accounting for stage_gain:
   total_DC_gain = total_DC_gain + AGC_params.AGC_stage_gain^(stage-1);
@@ -302,8 +303,10 @@
 AGC_coeffs.AGC_gain = total_DC_gain;
 
 % % print some results
-% AGC_coeffs
-% AGC_spatial_FIR = AGC_coeffs.AGC_spatial_FIR
+AGC_coeffs
+AGC_spatial_FIR = AGC_coeffs.AGC_spatial_FIR
+AGC_spatial_iterations = AGC_coeffs.AGC_spatial_iterations
+AGC_spatial_n_taps = AGC_coeffs.AGC_spatial_n_taps
 
 
 %%
@@ -453,7 +456,7 @@
 %                 AGC_polez2: [0.2219 0.3165 0.4260 0.4414]
 %     AGC_spatial_iterations: [1 1 2 2]
 %            AGC_spatial_FIR: [3x4 double]
-%                 AGC_n_taps: [3 5 5 5]
+%         AGC_spatial_n_taps: [3 5 5 5]
 %             AGC_mix_coeffs: [0 0.0454 0.0227 0.0113]
 %                   AGC_gain: 15
 %               detect_scale: 0.0664
--- a/trunk/matlab/bmm/carfac/CARFAC_Init.m	Fri Mar 16 04:31:56 2012 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_Init.m	Thu Mar 22 22:37:56 2012 +0000
@@ -28,7 +28,7 @@
 % level of object.  I like fewer structs and class types.
 
 if nargin < 2
-  n_ears = 1;  % monaural
+  n_ears = 1;  % monaural default
 end
 
 % % this is probably what I'd do in the C++ version:
@@ -54,38 +54,6 @@
   CF_struct.AGC_state(ear) = AGC_Init_State(CF_struct.AGC_coeffs);
 end
 
-% for ear = 1:n_ears
-%   CF_struct.CAR_state(ear).z1_memory = zeros(n_ch, 1);
-%   CF_struct.CAR_state(ear).z2_memory = zeros(n_ch, 1);
-%   CF_struct.CAR_state(ear).zA_memory = zeros(n_ch, 1);  % cubic loop
-%   CF_struct.CAR_state(ear).zB_memory = zeros(n_ch, 1);  % AGC interp
-%   CF_struct.CAR_state(ear).dzB_memory = zeros(n_ch, 1);  % AGC incr
-%   CF_struct.CAR_state(ear).zY_memory = zeros(n_ch, 1);
-%   CF_struct.CAR_state(ear).detect_accum = zeros(n_ch, 1);
-%   CF_struct.CAR_state(ear).g_memory = ...
-%     CF_struct.CAR_coeffs(ear).g0_coeffs;  % initial g for min_zeta
-%   CF_struct.CAR_state(ear).dg_memory = zeros(n_ch, 1);    % g interp
-%   % AGC loop filters' state:
-%   CF_struct.AGC_state(ear).AGC_memory = zeros(n_ch, n_AGC_stages);  % HACK init
-%   CF_struct.AGC_state(ear).input_accum = zeros(n_ch, n_AGC_stages);  % HACK init
-%   % IHC state:
-%   if CF_struct.IHC_coeffs.just_hwr
-%     CF_struct.IHC_state(ear).ihc_accum = zeros(n_ch, 1);
-%   else
-%     CF_struct.IHC_state(ear).cap_voltage = ...
-%       CF_struct.IHC_coeffs.rest_cap * ones(n_ch, 1);
-%     CF_struct.IHC_state(ear).cap1_voltage = ...
-%       CF_struct.IHC_coeffs.rest_cap1 * ones(n_ch, 1);
-%     CF_struct.IHC_state(ear).cap2_voltage = ...
-%       CF_struct.IHC_coeffs.rest_cap2 * ones(n_ch, 1);
-%     CF_struct.IHC_state(ear).lpf1_state = ...
-%       CF_struct.IHC_coeffs.rest_output * zeros(n_ch, 1);
-%     CF_struct.IHC_state(ear).lpf2_state = ...
-%       CF_struct.IHC_coeffs.rest_output * zeros(n_ch, 1);
-%     CF_struct.IHC_state(ear).ihc_accum = zeros(n_ch, 1);
-%   end
-% end
-
 
 function state = CAR_Init_State(coeffs)
 n_ch = coeffs.n_ch;
@@ -96,7 +64,6 @@
   'zB_memory', zeros(n_ch, 1), ...
   'dzB_memory', zeros(n_ch, 1), ...
   'zY_memory', zeros(n_ch, 1), ...
-  'detect_accum', zeros(n_ch, 1), ...
   'g_memory', coeffs.g0_coeffs, ...
   'dg_memory', zeros(n_ch, 1) ...
   );
--- a/trunk/matlab/bmm/carfac/CARFAC_Run.m	Fri Mar 16 04:31:56 2012 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_Run.m	Thu Mar 22 22:37:56 2012 +0000
@@ -17,7 +17,7 @@
 % See the License for the specific language governing permissions and
 % limitations under the License.
 
-function [CF, decim_naps, naps] = CARFAC_Run ...
+function [CF, decim_naps, naps, BM] = CARFAC_Run ...
   (CF, input_waves, AGC_plot_fig_num)
 % function [CF, decim_naps, naps] = CARFAC_Run ...
 %   (CF, input_waves, AGC_plot_fig_num)
@@ -50,6 +50,12 @@
   AGC_plot_fig_num = 0;
 end
 
+if nargout > 3
+  BM = zeros(n_samp, n_ch, n_ears);
+else
+  BM = [];
+end
+
 if n_ears ~= CF.n_ears
   error('bad number of input_waves channels passed to CARFAC_Run')
 end
@@ -57,7 +63,7 @@
 
 naps = zeros(n_samp, n_ch, n_ears);
 
-seglen = 16;
+seglen = 256;
 n_segs = ceil(n_samp / seglen);
 
 if nargout > 1
@@ -82,7 +88,18 @@
     k_range = seglen*(seg_num - 1) + (1:seglen);
   end
   % Process a segment to get a slice of decim_naps, and plot AGC state:
-  [seg_naps, CF] = CARFAC_Run_Segment(CF, input_waves(k_range, :));
+  if ~isempty(BM)
+    [seg_naps, CF, seg_BM] = CARFAC_Run_Segment(CF, input_waves(k_range, :));
+  else
+    [seg_naps, CF] = CARFAC_Run_Segment(CF, input_waves(k_range, :));
+  end
+  
+  if ~isempty(BM)
+    for ear = 1:n_ears
+      % Accumulate segment BM to make full BM
+      BM(k_range, :, ear) = seg_BM(:, :, ear);
+    end
+  end
   
   if ~isempty(naps)
     for ear = 1:n_ears
@@ -109,11 +126,14 @@
       for stage = 1:3;
         plot(2^(stage-1) * (CF.AGC_state(ear).AGC_memory(:, stage) - ...
           2 * CF.AGC_state(ear).AGC_memory(:, stage+1)));
+        plot(2^(stage-1) * (CF.AGC_state(ear).AGC_memory(:, stage) - ...
+          0 * CF.AGC_state(ear).AGC_memory(:, stage+1)), 'r--');
       end
       stage = 4;
       plot(2^(stage-1) * CF.AGC_state(ear).AGC_memory(:, stage));
+      plot(2^(stage-1) * CF.AGC_state(ear).AGC_memory(:, stage), 'r--');
     end
-    axis([0, CF.n_ch+1, 0.0, max(maxes) + 0.01]);
+    axis([0, CF.n_ch+1, 0.0, max(maxes) * 1.01 + 0.002]);
     drawnow
   end
 
--- a/trunk/matlab/bmm/carfac/CARFAC_Run_Linear.m	Fri Mar 16 04:31:56 2012 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_Run_Linear.m	Thu Mar 22 22:37:56 2012 +0000
@@ -49,7 +49,7 @@
 for k = 1:n_samp
   % at each time step, possibly handle multiple channels
   for ear = 1:n_ears
-    [filters_out, CF.CAR_state(ear)] = CARFAC_FilterStep( ...
+    [filters_out, CF.CAR_state(ear)] = CARFAC_CAR_Step( ...
       input_waves(k, ear), CF.CAR_coeffs, CF.CAR_state(ear));
     naps(k, :, ear) = filters_out;  % linear
   end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_Run_Open_Loop.m	Thu Mar 22 22:37:56 2012 +0000
@@ -0,0 +1,119 @@
+% Copyright 2012, Google, Inc.
+% Author Richard F. Lyon
+%
+% This Matlab file is part of an implementation of Lyon's cochlear model:
+% "Cascade of Asymmetric Resonators with Fast-Acting Compression"
+% to supplement Lyon's upcoming book "Human and Machine Hearing"
+%
+% Licensed under the Apache License, Version 2.0 (the "License");
+% you may not use this file except in compliance with the License.
+% You may obtain a copy of the License at
+%
+%     http://www.apache.org/licenses/LICENSE-2.0
+%
+% Unless required by applicable law or agreed to in writing, software
+% distributed under the License is distributed on an "AS IS" BASIS,
+% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+% See the License for the specific language governing permissions and
+% limitations under the License.
+
+function [CF, decim_naps, naps] = CARFAC_Run_Open_Loop ...
+  (CF, input_waves, AGC_plot_fig_num)
+% function [CF, decim_naps, naps] = CARFAC_Run_Open_Loop ...
+%   (CF, input_waves, AGC_plot_fig_num)
+%
+% Freeze the damping by disabling AGC feedback, and run so we can
+% see what the filters and AGC do in that frozen state.  And zap the
+% stage gain in the AGC so we can see the state filters without combining
+% them.
+
+[n_samp, n_ears] = size(input_waves);
+n_ch = CF.n_ch;
+
+if nargin < 3
+  AGC_plot_fig_num = 0;
+end
+
+if n_ears ~= CF.n_ears
+  error('bad number of input_waves channels passed to CARFAC_Run')
+end
+
+
+naps = zeros(n_samp, n_ch, n_ears);
+
+seglen = 16;
+n_segs = ceil(n_samp / seglen);
+
+if nargout > 1
+  % make decimated detect output:
+  decim_naps = zeros(n_segs, CF.n_ch, CF.n_ears);
+else
+  decim_naps = [];
+end
+
+if nargout > 2
+  % make decimated detect output:
+  naps = zeros(n_samp, CF.n_ch, CF.n_ears);
+else
+  naps = [];
+end
+
+% zero the deltas:
+for ear = 1:CF.n_ears
+  CF.CAR_state(ear).dzB_memory = 0;
+  CF.CAR_state(ear).dg_memory = 0;
+end
+open_loop = 1;
+CF.AGC_coeffs.AGC_stage_gain = 0;  % HACK to see the stages separately
+
+smoothed_state = 0;
+
+for seg_num = 1:n_segs
+  if seg_num == n_segs
+    % The last segement may be short of seglen, but do it anyway:
+    k_range = (seglen*(seg_num - 1) + 1):n_samp;
+  else
+    k_range = seglen*(seg_num - 1) + (1:seglen);
+  end
+  % Process a segment to get a slice of decim_naps, and plot AGC state:
+  [seg_naps, CF] = CARFAC_Run_Segment(CF, input_waves(k_range, :), ...
+    open_loop);
+  
+  if ~isempty(naps)
+    for ear = 1:n_ears
+      % Accumulate segment naps to make full naps
+      naps(k_range, :, ear) = seg_naps(:, :, ear);
+    end
+  end
+  
+  if ~isempty(decim_naps)
+    for ear = 1:n_ears
+      decim_naps(seg_num, :, ear) = CF.IHC_state(ear).ihc_accum / seglen;
+      CF.IHC_state(ear).ihc_accum = zeros(n_ch,1);
+    end
+  end
+  
+  if AGC_plot_fig_num
+    figure(AGC_plot_fig_num); hold off; clf
+    set(gca, 'Position', [.25, .25, .5, .5])
+    smoothed_state = (3*smoothed_state + CF.AGC_state(1).AGC_memory) / 4;
+    for ear = 1
+      total_state = 0;
+      for stage = 1:4;
+        weighted_state = smoothed_state(:, stage) * 2^(stage-1);
+        plot(weighted_state, 'k-', 'LineWidth', 0.4);
+        hold on
+        total_state = total_state + weighted_state;
+      end
+      maxes(ear) = max(total_state);
+      plot(total_state, 'k-', 'LineWidth', 1.1)
+    end
+    
+    axis([0, CF.n_ch+1, 0.0, max(maxes) * 1.01 + 0.002]);
+    drawnow
+  end
+  
+end
+
+
+
--- a/trunk/matlab/bmm/carfac/CARFAC_Run_Segment.m	Fri Mar 16 04:31:56 2012 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_Run_Segment.m	Thu Mar 22 22:37:56 2012 +0000
@@ -17,8 +17,8 @@
 % See the License for the specific language governing permissions and
 % limitations under the License.
 
-function [naps, CF] = CARFAC_Run_Segment(CF, input_waves)
-% function [naps, CF, decim_naps] = CARFAC_Run_Segment(CF, input_waves)
+function [naps, CF, BM] = CARFAC_Run_Segment(CF, input_waves, open_loop)
+% function [naps, CF, BM] = CARFAC_Run_Segment(CF, input_waves, open_loop)
 % 
 % This function runs the CARFAC; that is, filters a 1 or more channel
 % sound input segment to make one or more neural activity patterns (naps);
@@ -31,7 +31,7 @@
 %
 % naps has a row per time sample, a column per filterbank channel, and
 % a layer per audio channel if more than 1.
-% decim_naps is like naps but time-decimated by the int CF.decimation.
+% BM is basilar membrane motion (filter outputs before detection).
 %
 % the input_waves are assumed to be sampled at the same rate as the
 % CARFAC is designed for; a resampling may be needed before calling this.
@@ -45,6 +45,16 @@
 % CF = CARFAC_Design(fs, CF_CAR_params, CF_AGC_params, n_ears)
 % transfns = CARFAC_Transfer_Functions(CF, to_chans, from_chans)
 
+if nargin < 3
+  open_loop = 0;
+end
+
+if nargout > 2
+  do_BM = 1;
+else
+  do_BM = 0;
+end
+
 [n_samp, n_ears] = size(input_waves);
 
 if n_ears ~= CF.n_ears
@@ -53,6 +63,9 @@
 
 n_ch = CF.n_ch;
 naps = zeros(n_samp, n_ch, n_ears);  % allocate space for result
+if do_BM
+  BM = zeros(n_samp, n_ch, n_ears);
+end
 
 detects = zeros(n_ch, n_ears);
 for k = 1:n_samp
@@ -67,6 +80,9 @@
     
     detects(:, ear) = ihc_out;  % for input to AGC, and out to SAI
     naps(k, :, ear) = ihc_out;  % output to neural activity pattern  
+    if do_BM
+      BM(k, :, ear) = car_out;
+    end
   end
   % run the AGC update step, taking input from IHC_state, 
   % decimating internally, all ears at once due to mixing across them:
@@ -74,7 +90,7 @@
     CF.AGC_coeffs, detects, CF.AGC_state);
   
   % connect the feedback from AGC_state to CAR_state when it updates
-  if updated
+  if updated & ~open_loop
     CF = CARFAC_Close_AGC_Loop(CF);
   end
 end
--- a/trunk/matlab/bmm/carfac/CARFAC_Spatial_Smooth.m	Fri Mar 16 04:31:56 2012 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_Spatial_Smooth.m	Thu Mar 22 22:37:56 2012 +0000
@@ -27,7 +27,7 @@
 
 if use_FIR
   FIR_coeffs = coeffs.AGC_spatial_FIR(:,stage);
-  switch coeffs.AGC_n_taps(stage)
+  switch coeffs.AGC_spatial_n_taps(stage)
     case 3
       for iter = 1:n_iterations
         stage_state = ...
@@ -45,10 +45,10 @@
           stage_state([3:end, end, end-1], :));
       end
     otherwise
-      error('Bad n_taps in CARFAC_Spatial_Smooth');
+      error('Bad AGC_spatial_n_taps in CARFAC_Spatial_Smooth');
   end
 else
-  % use IIR method, back-and-forth firt-order smoothers:
+  % use IIR method, back-and-forth first-order smoothers:
   stage_state = SmoothDoubleExponential(stage_state, ...
     coeffs.AGC_polez1(stage), coeffs.AGC_polez2(stage));
 end
--- a/trunk/matlab/bmm/carfac/CARFAC_hacking.m	Fri Mar 16 04:31:56 2012 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_hacking.m	Thu Mar 22 22:37:56 2012 +0000
@@ -22,17 +22,33 @@
 clear variables
 
 %%
-file_signal = wavread('plan.wav');
-
-% file_signal = file_signal(9000+(1:10000));  % trim for a faster test
-file_signal = file_signal(9300+(1:5000));  % trim for a faster test
+use_plan_file = 0;
+if use_plan_file
+  
+  file_signal = wavread('plan.wav');
+  file_signal = file_signal(8100+(1:20000));  % trim for a faster test
+  
+else
+    flist = [1000];
+    alist = [1];
+    flist = 1000;
+    alist = 1;
+    sine_signal = 0;
+    times = (0:19999)' / 22050;
+    for fno = 1:length(flist)
+      sine_signal = sine_signal + alist(fno)*sin(flist(fno)*2*pi*times);
+    end
+    growth_power = 0;  % use 0 for flat, 4 or more for near exponential
+    file_signal = 1.0 * (sine_signal .* (times/max(times)).^growth_power);
+end
 
 % repeat with negated signal to compare responses:
 % file_signal = [file_signal; -file_signal];
 
 % make a long test signal by repeating at different levels:
-test_signal = [];
-for dB = -40:20:0  % -60:20:40  % -80:20:60
+dB = -80;
+test_signal =  10^(dB/20)* file_signal(1:4000) % lead-in [];
+for dB =  -80:20:80
   test_signal = [test_signal; file_signal * 10^(dB/20)];
 end
 
@@ -43,16 +59,20 @@
 
 agc_plot_fig_num = 6;
 
-for n_ears = 1:2
+for n_ears = 1  % 1:2
   CF_struct = CARFAC_Init(CF_struct, n_ears);
 
-  [CF_struct, nap_decim, nap] = CARFAC_Run(CF_struct, test_signal, ...
+  [CF_struct, nap_decim, nap, BM] = CARFAC_Run(CF_struct, test_signal, ...
     agc_plot_fig_num);
 
 %   nap = deskew(nap);  % deskew doesn't make much difference
 
+%   dB_BM = 10/log(10) * log(filter(1, [1, -0.995], BM(:, 38:40, :).^2));
+  dB_BM = 10/log(10) * log(filter(1, [1, -0.995], BM(:, 20:50, :).^2));
+
   if n_ears == 1  % because this hack doesn't work for binarual yet
-    MultiScaleSmooth(nap_decim, 4);
+    MultiScaleSmooth(dB_BM(5000:200:end, :, :), 1);
+%     MultiScaleSmooth(nap_decim, 4);
   end
 
   % Display results for 1 or 2 ears:
--- a/trunk/matlab/bmm/carfac/MultiScaleSmooth.m	Fri Mar 16 04:31:56 2012 +0000
+++ b/trunk/matlab/bmm/carfac/MultiScaleSmooth.m	Thu Mar 22 22:37:56 2012 +0000
@@ -61,10 +61,3 @@
 
 end
 
-
-function waves = deskew(waves)
-
-for col = 1:size(waves, 2)
-  waves(1:(end-col+1), col) = waves(col:end, col);
-end
-