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