dicklyon@671: % Copyright 2010, Google, Inc.
dicklyon@671: % Author Richard F. Lyon
dicklyon@671: %
dicklyon@671: % This Matlab file is part of an implementation of Lyon's cochlear model:
dicklyon@671: % "Cascade of Asymmetric Resonators with Fast-Acting Compression"
dicklyon@671: % to supplement Lyon's upcoming book "Human and Machine Hearing"
dicklyon@671: %
dicklyon@671: % Licensed under the Apache License, Version 2.0 (the "License");
dicklyon@671: % you may not use this file except in compliance with the License.
dicklyon@671: % You may obtain a copy of the License at
dicklyon@671: %
dicklyon@671: %     http://www.apache.org/licenses/LICENSE-2.0
dicklyon@671: %
dicklyon@671: % Unless required by applicable law or agreed to in writing, software
dicklyon@671: % distributed under the License is distributed on an "AS IS" BASIS,
dicklyon@671: % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
dicklyon@671: % See the License for the specific language governing permissions and
dicklyon@671: % limitations under the License.
dicklyon@671: 
dicklyon@671: function cols = smooth1d(cols, scale)
dicklyon@671: % smooth1d - Smooth the columns of the input.
dicklyon@671: %
dicklyon@671: % output = smooth1d(input, smoothing_factor)
dicklyon@671: %
dicklyon@671: % Smooth the columns of input using a one-pole smoothing filter, using the
dicklyon@671: % provided smoothing factor.
dicklyon@671: %
dicklyon@671: % TODO(dross, dicklyon): make this code satisfy the google3 Matlab style.
dicklyon@671: 
dicklyon@671: [nr, nc, nl] = size(cols);
dicklyon@671: if nr == 1
dicklyon@671:     if nl == 1
dicklyon@671:         cols = cols';
dicklyon@671:         [nr, nc, nl] = size(cols);
dicklyon@671:     else
dicklyon@671:         disp('error in shape passed to smooth1d')
dicklyon@671:     end
dicklyon@671: end
dicklyon@671: 
dicklyon@671: if scale==0
dicklyon@671:     polez = 0; % no smoothing at all
dicklyon@671:     return;
dicklyon@671: else
dicklyon@671:     % this coefficient matches the curvature at DC of a discrete Gaussian:
dicklyon@671:     t = scale^2;
dicklyon@671:     polez = 1 + 1/t - sqrt((1+1/t)^2 - 1); % equiv. to the one below
dicklyon@671:     % polez = 1 + 1/t - sqrt((2/t) + 1/t^2); % equiv. to the one above
dicklyon@671:     % polez is Z position of real pole
dicklyon@671: end
dicklyon@671: 
dicklyon@671: [x, state] = filter(1-polez, [1, -polez], cols);
dicklyon@671: cols = cols(end:-1:1,:);
dicklyon@671: [cols, state] = filter(1-polez, [1, -polez], cols, state);
dicklyon@671: cols = cols(end:-1:1,:);
dicklyon@671: [cols, state] = filter(1-polez, [1, -polez], cols, state);
dicklyon@671: 
dicklyon@671: return