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 signal_vecs = SmoothDoubleExponential(signal_vecs, ...
tom@516:   polez1, polez2, fast_matlab_way)
tom@516: % function signal_vecs = SmoothDoubleExponential(signal_vecs, ...
tom@516: %   polez1, polez2, fast_matlab_way)
tom@516: %
tom@516: % Smooth the input column vectors in signal_vecs using forward
tom@516: % and backwards one-pole smoothing filters, backwards first, with
tom@516: % approximately reflecting edge conditions.
tom@516: %
tom@516: % It will be done with Matlab's filter function if "fast_matlab_way"
tom@516: % is nonzero or defaulted; use 0 to test the algorithm for how to do it
tom@516: % in sequential c code.
tom@516: 
tom@516: if nargin < 4
tom@516:   fast_matlab_way = 1;
tom@516:   % can also use the slow way with explicit loop like we'll do in C++
tom@516: end
tom@516: 
tom@516: if fast_matlab_way
tom@516:   [junk, Z_state] = filter(1-polez1, [1, -polez1], ...
tom@516:     signal_vecs((end-10):end, :));  % initialize state from 10 points
tom@516:   [signal_vecs(end:-1:1), Z_state] = filter(1-polez2, [1, -polez2], ...
tom@516:     signal_vecs(end:-1:1), Z_state*polez2/polez1);
tom@516:   signal_vecs = filter(1-polez1, [1, -polez1], signal_vecs, ...
tom@516:     Z_state*polez1/polez2);
tom@516: else
tom@516:   npts = size(signal_vecs, 1);
tom@516:   state = zeros(size(signal_vecs, 2));
tom@516:   for index = npts-10:npts
tom@516:     input = signal_vecs(index, :);
tom@516:     state = state + (1 - polez1) * (input - state);
tom@516:   end
tom@516:   % smooth backward with polez2, starting with state from above:
tom@516:   for index = npts:-1:1
tom@516:     input = signal_vecs(index, :);
tom@516:     state = state + (1 - polez2) * (input - state);
tom@516:     signal_vecs(index, :) = state;
tom@516:   end
tom@516:   % smooth forward with polez1, starting with state from above:
tom@516:   for index = 1:npts
tom@516:     input = signal_vecs(index, :);
tom@516:     state = state + (1 - polez1) * (input - state);
tom@516:     signal_vecs(index, :) = state;
tom@516:   end
tom@516: end
dicklyon@523: