Mercurial > hg > aimc

--- a/trunk/matlab/bmm/carfac/CARFAC_CAR_Step.m	Tue May 01 23:32:24 2012 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_CAR_Step.m	Sat May 12 04:31:59 2012 +0000
@@ -29,9 +29,6 @@
 zB = state.zB_memory + state.dzB_memory; % AGC interpolation
 r1 = CAR_coeffs.r1_coeffs;
 g = state.g_memory + state.dg_memory;  % interp g
-v_offset  = CAR_coeffs.v_offset;
-v2_corner = CAR_coeffs.v2_corner;
-v_damp_max = CAR_coeffs.v_damp_max;

 % zB and zA are "extra damping", and multiply zr (compressed theta):
 r = r1 - CAR_coeffs.zr_coeffs .* (zA + zB);
@@ -43,11 +40,8 @@
 z2 = r .* (CAR_coeffs.c0_coeffs .* state.z1_memory + ...
   CAR_coeffs.a0_coeffs .* state.z2_memory);

-% update the "velocity" for cubic nonlinearity, into zA:
-zA = (((state.z2_memory - z2) .* CAR_coeffs.velocity_scale) + ...
-  v_offset) .^ 2;
-% soft saturation to make it more like an "essential" nonlinearity:
-zA = v_damp_max * zA ./ (v2_corner + zA);
+% update the nonlinear function of "velocity", into zA:
+zA = CARFAC_OHC_NLF(state.z2_memory - z2, CAR_coeffs);

 zY = CAR_coeffs.h_coeffs .* z2;  % partial output
--- a/trunk/matlab/bmm/carfac/CARFAC_Design.m	Tue May 01 23:32:24 2012 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_Design.m	Sat May 12 04:31:59 2012 +0000
@@ -53,7 +53,7 @@
 if nargin < 3
   CF_CAR_params = struct( ...
     'velocity_scale', 0.2, ...  % for the "cubic" velocity nonlinearity
-    'v_offset', 0.01, ...  % offset gives a quadratic part
+    'v_offset', 0.04, ...  % offset gives a quadratic part
     'v2_corner', 0.2, ...  % corner for essential nonlin
     'v_damp_max', 0.01, ... % damping delta damping from velocity nonlin
     'min_zeta', 0.10, ... % minimum damping factor in mid-freq channels
@@ -71,7 +71,7 @@
     'n_stages', 4, ...
     '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
+    'decimation', [4, 2, 2, 2], ...  % how often to update the AGC states
     '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
--- a/trunk/matlab/bmm/carfac/CARFAC_Design_Doc.txt	Tue May 01 23:32:24 2012 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_Design_Doc.txt	Sat May 12 04:31:59 2012 +0000
@@ -7,7 +7,7 @@
 processor, for mono, stereo, or multi-channel sound inputs.  This file
 describes aspects of the software design.

-The implementation will include equivalent Maltab, C++, Python, and
+The implementation will include equivalent Matlab, C++, Python, and
 perhaps other versions.  They should all be based on the same set of
 classes (or structs) and functions, with roughly equivalent
 functionality and names where possible.  The examples here are Matlab,
@@ -20,7 +20,10 @@
 includes various sub-objects representing the parameters, designs, and
 states of different parts of the CAR-FAC model.

-The three main subsystems of the CARFAC are the cascade of asymmetric
+The CARFAC class includes a vector of "ear" objects -- one for mono,
+two for stereo, or more.
+
+The three main subsystems of the EAR are the cascade of asymmetric
 resonators (CAR), the inner hair cell (IHC), and the automatic gain
 control (AGC).  These are not intended to work independently, so may
 not have corresponding classes to combine all their aspects (in the
@@ -35,8 +38,6 @@
 These names can be used both for the classes, and slightly modified
 for the member variables, arguments, temps, etc.

- -- (presently we have "filters" where I say "CAR"; I'll change it) --
-
 The params are inputs that specify the system; the coeffs are things
 like filter coefficients, things computed once and used at run time;
 the state is whatever internal state is needed between running the
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_OHC_NLF.m	Sat May 12 04:31:59 2012 +0000
@@ -0,0 +1,7 @@
+function nlf = CARFAC_OHC_NLF(velocities, CAR_coeffs)
+
+nlf = ((velocities .* CAR_coeffs.velocity_scale) + ...
+  CAR_coeffs.v_offset) .^ 2;
+% soft saturation to make it more like an "essential" nonlinearity:
+nlf = CAR_coeffs.v_damp_max * nlf ./ (CAR_coeffs.v2_corner + nlf);
+
--- a/trunk/matlab/bmm/carfac/CARFAC_Run.m	Tue May 01 23:32:24 2012 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_Run.m	Sat May 12 04:31:59 2012 +0000
@@ -17,9 +17,9 @@
 % See the License for the specific language governing permissions and
 % limitations under the License.

-function [CF, decim_naps, naps, BM] = CARFAC_Run ...
+function [CF, decim_naps, naps, BM, ohc, agc] = CARFAC_Run ...
   (CF, input_waves, AGC_plot_fig_num)
-% function [CF, decim_naps, naps] = CARFAC_Run ...
+% function [CF, decim_naps, naps, BM, ohc, agc] = CARFAC_Run ...
 %   (CF, input_waves, AGC_plot_fig_num)
 % This function runs the CARFAC; that is, filters a 1 or more channel
 % sound input to make one or more neural activity patterns (naps).
@@ -38,10 +38,7 @@
 % 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.
 %
-% The function works as an outer iteration on time, updating all the
-% filters and AGC states concurrently, so that the different channels can
-% interact easily.  The inner loops are over filterbank channels, and
-% this level should be kept efficient.
+% ohc and agc are optional extra outputs for diagnosing internals.

 [n_samp, n_ears] = size(input_waves);
 n_ch = CF.n_ch;
@@ -56,6 +53,18 @@
   BM = [];
 end

+if nargout > 4
+  ohc = zeros(n_samp, n_ch, n_ears);
+else
+  ohc = [];
+end
+
+if nargout > 5
+  agc = zeros(n_samp, n_ch, n_ears);
+else
+  agc = [];
+end
+
 if n_ears ~= CF.n_ears
   error('bad number of input_waves channels passed to CARFAC_Run')
 end
@@ -89,7 +98,8 @@
   end
   % Process a segment to get a slice of decim_naps, and plot AGC state:
   if ~isempty(BM)
-    [seg_naps, CF, seg_BM] = CARFAC_Run_Segment(CF, input_waves(k_range, :));
+    % ask for everything in this case, for laziness:
+    [seg_naps, CF, seg_BM, seg_ohc, seg_agc] = CARFAC_Run_Segment(CF, input_waves(k_range, :));
   else
     [seg_naps, CF] = CARFAC_Run_Segment(CF, input_waves(k_range, :));
   end
@@ -108,6 +118,20 @@
     end
   end

+  if ~isempty(ohc)
+    for ear = 1:n_ears
+      % Accumulate segment naps to make full naps
+      ohc(k_range, :, ear) = seg_ohc(:, :, ear);
+    end
+  end
+
+  if ~isempty(agc)
+    for ear = 1:n_ears
+      % Accumulate segment naps to make full naps
+      agc(k_range, :, ear) = seg_agc(:, :, ear);
+    end
+  end
+
   if ~isempty(decim_naps)
     for ear = 1:n_ears
       decim_naps(seg_num, :, ear) = CF.ears(ear).IHC_state.ihc_accum / seglen;
--- a/trunk/matlab/bmm/carfac/CARFAC_Run_Segment.m	Tue May 01 23:32:24 2012 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_Run_Segment.m	Sat May 12 04:31:59 2012 +0000
@@ -17,8 +17,10 @@
 % See the License for the specific language governing permissions and
 % limitations under the License.

-function [naps, CF, BM] = CARFAC_Run_Segment(CF, input_waves, open_loop)
-% function [naps, CF, BM] = CARFAC_Run_Segment(CF, input_waves, open_loop)
+function [naps, CF, BM, seg_ohc, seg_agc] = CARFAC_Run_Segment(...
+  CF, input_waves, open_loop)
+% function [naps, CF, BM, seg_ohc, seg_agc] = 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);
@@ -41,9 +43,8 @@
 % interact easily.  The inner loops are over filterbank channels, and
 % this level should be kept efficient.
 %
-% See other functions for designing and characterizing the CARFAC:
-% CF = CARFAC_Design(fs, CF_CAR_params, CF_AGC_params, n_ears)
-% transfns = CARFAC_Transfer_Functions(CF, to_chans, from_chans)
+% seg_ohc seg_agc are optional extra outputs useful for seeing what the
+% ohc nonlinearity and agc are doing; both in terms of extra damping.

 if nargin < 3
   open_loop = 0;
@@ -65,12 +66,18 @@
 naps = zeros(n_samp, n_ch, n_ears);  % allocate space for result
 if do_BM
   BM = zeros(n_samp, n_ch, n_ears);
+  seg_ohc = zeros(n_samp, n_ch, n_ears);
+  seg_agc = zeros(n_samp, n_ch, n_ears);
 end

 detects = zeros(n_ch, n_ears);
 for k = 1:n_samp
   % at each time step, possibly handle multiple channels
   for ear = 1:n_ears
+
+    % This would be cleaner if we could just get and use a reference to
+    % CF.ears(ear), but Matlab doesn't work that way...
+
     [car_out, CF.ears(ear).CAR_state] = CARFAC_CAR_Step( ...
       input_waves(k, ear), CF.ears(ear).CAR_coeffs, CF.ears(ear).CAR_state);

@@ -86,6 +93,9 @@
     naps(k, :, ear) = ihc_out;  % output to neural activity pattern
     if do_BM
       BM(k, :, ear) = car_out;
+      state = CF.ears(ear).CAR_state;
+      seg_ohc(k, :, ear) = state.zA_memory;
+      seg_agc(k, :, ear) = state.zB_memory;;
     end
   end
--- a/trunk/matlab/bmm/carfac/CARFAC_hacking.m	Tue May 01 23:32:24 2012 +0000
+++ b/trunk/matlab/bmm/carfac/CARFAC_hacking.m	Sat May 12 04:31:59 2012 +0000
@@ -26,29 +26,37 @@
 if use_plan_file

   file_signal = wavread('plan.wav');
-  file_signal = file_signal(8100+(1:20000));  % trim for a faster test
+  %   file_signal = file_signal(8100+(1:20000));  % trim for a faster test
+  file_signal = file_signal(10000+(1:10000));  % 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);
+  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:
-dB = -80;
-test_signal =  10^(dB/20)* file_signal(1:4000) % lead-in [];
-for dB =  -80:20:60
+do_distortion_figs = 0;  % use 1 for distortion figure hack....
+if do_distortion_figs
+  dB_list = -40:20:0
+else
+  dB_list = -80:20:60
+end
+% dB = dB_list(1);
+% test_signal =  10^(dB/20)* file_signal(1:4000) % lead-in [];
+test_signal = [];
+for dB =  dB_list
   test_signal = [test_signal; file_signal * 10^(dB/20)];
 end

@@ -57,28 +65,28 @@

 agc_plot_fig_num = 6;

-for n_ears = 1:2
+for n_ears = 1    %  1:2

   CF_struct = CARFAC_Design(n_ears);  % default design
-
+
   if n_ears == 2
     % For the 2-channel pass, add a silent second channel:
     test_signal = [test_signal, zeros(size(test_signal))];
   end

   CF_struct = CARFAC_Init(CF_struct);
-
-  [CF_struct, nap_decim, nap, BM] = CARFAC_Run(CF_struct, test_signal, ...
+
+  [CF_struct, nap_decim, nap, BM, ohc, agc] = 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));
-
+
+  %   nap = deskew(nap);  % deskew doesn't make much difference
+
+  %   dB_BM = 10/log(10) * log(filter(1, [1, -0.995], BM(:, 20:50, :).^2));
+  dB_BM = 10/log(10) * log(filter(1, [1, -0.995], BM(:, :, :).^2));
+
   % only ear 1:
   MultiScaleSmooth(dB_BM(5000:200:end, :, 1), 1);
-
+
   % Display results for 1 or 2 ears:
   for ear = 1:n_ears
     smooth_nap = nap_decim(:, :, ear);
@@ -90,13 +98,59 @@
     title('smooth nap from nap decim')
     colormap(1 - gray);
   end
-
+
   % Show resulting data, even though M-Lint complains:
   CF_struct
-  CF_struct.CAR_state
-  CF_struct.AGC_state
+  CF_struct.ears(1).CAR_state
+  CF_struct.ears(1).AGC_state
   min_max_decim = [min(nap_decim(:)), max(nap_decim(:))]
-
+
+  %%
+  if do_distortion_figs
+    channels = [38, 39, 40];
+    times = 0000 + (3200:3599);
+    smoothed_ohc = ohc;
+    epsi = 0.4;
+    %   smoothed_ohc = filter(epsi, [1, -(1-epsi)], smoothed_ohc);
+    %   smoothed_ohc = filter(epsi, [1, -(1-epsi)], smoothed_ohc);
+
+    figure(101); hold off
+    plot(smoothed_ohc(times, channels))
+    hold on
+    plot(smoothed_ohc(times, channels(2)), 'k-', 'LineWidth', 1.5)
+    title('OHC')
+
+    figure(105); hold off
+    plot(agc(times, channels))
+    hold on
+    plot(agc(times, channels(2)), 'k-', 'LineWidth', 1.5)
+    title('AGC')
+
+    figure(103); hold off
+    plot(BM(times, channels))
+    hold on
+    plot(BM(times, channels(2)), 'k-', 'LineWidth', 1.5)
+    title('BM')
+
+    extra_damping = smoothed_ohc + agc;
+    figure(102); hold off
+    plot(extra_damping(times, channels))
+    hold on
+    plot(extra_damping(times, channels(2)), 'k-', 'LineWidth', 1.5)
+    title('extra damping')
+
+    distortion = -(extra_damping - smooth1d(extra_damping, 10)) .* BM;
+    distortion = filter(epsi, [1, -(1-epsi)], distortion);
+    distortion = filter(epsi, [1, -(1-epsi)], distortion);
+    %   distortion = filter(epsi, [1, -(1-epsi)], distortion);
+    figure(104); hold off
+    plot(distortion(times, channels))
+    hold on
+    plot(distortion(times, channels(2)), 'k-', 'LineWidth', 1.5)
+    title('distortion')
+
+  end
+  %%
 end

 % Expected result:  Figure 3 looks like figure 2, a tiny bit darker.