Daniel@0: function y=LyonPassiveEar(x,sr,df,earQ,stepfactor,differ,agcf,taufactor)
Daniel@0: % y=LyonPassiveEar(input, sample_rate, decimation, earQ, stepfactor);
Daniel@0: %
Daniel@0: % Changes
Daniel@0: % Zeroed out the preemphasis channels before doing the AGC.  This is needed
Daniel@0: % so that cochlear inversion is possible (since we don't have those channels,
Daniel@0: % we can't invert their gain.)..... Malcolm 6/26/95
Daniel@0: 
Daniel@0: % (c) 1998 Interval Research Corporation  
Daniel@0: 
Daniel@0: if nargin < 3
Daniel@0: 	fprintf('Syntax: y=LyonPassiveEar(input, sample_rate, decimation, earQ, stepfactor)\n');
Daniel@0: 	fprintf(' The input, sample_rate, and decimation parameters are mandatory.\n');
Daniel@0: 	return;
Daniel@0: end
Daniel@0: 
Daniel@0: if df < 1; df = 1; end
Daniel@0: if nargin < 4; earQ = 8; end
Daniel@0: if nargin < 5; stepfactor = earQ/32; end
Daniel@0: if nargin < 6; differ=1; end
Daniel@0: if nargin < 7; agcf=1; end
Daniel@0: if nargin < 8; taufactor=3; end
Daniel@0: 
Daniel@0: earFilters = DesignLyonFilters(sr, earQ, stepfactor);
Daniel@0: 
Daniel@0: nSamples = length(x);
Daniel@0: nOutputSamples = floor(nSamples/df);
Daniel@0: [nChannels filterWidth] = size(earFilters);
Daniel@0: 
Daniel@0: sosOutput = zeros(nChannels, df);
Daniel@0: sosState = zeros(nChannels, 2);
Daniel@0: agcState = zeros(nChannels, 4);
Daniel@0: y = zeros(nChannels, nOutputSamples);
Daniel@0: 
Daniel@0: decEps = EpsilonFromTauFS(df/sr*taufactor,sr);
Daniel@0: decState = zeros(nChannels, 2);
Daniel@0: decFilt = SetGain([0 0 1 -2*(1-decEps) (1-decEps)^2], 1, 0, sr);
Daniel@0: 
Daniel@0: eps1 = EpsilonFromTauFS(.64,sr);
Daniel@0: eps2 = EpsilonFromTauFS(.16,sr);
Daniel@0: eps3 = EpsilonFromTauFS(.04,sr);
Daniel@0: eps4 = EpsilonFromTauFS(.01,sr);
Daniel@0: 
Daniel@0: tar1 = .0032;
Daniel@0: tar2 = .0016;
Daniel@0: tar3 = .0008;
Daniel@0: tar4 = .0004;
Daniel@0: 
Daniel@0: if 0
Daniel@0: 	fprintf('df=%g, earq=%g, stepfactor=%g, differ=%g\n',df,earQ,stepfactor,differ);
Daniel@0: 	fprintf('agcf=%g, taufactor=%g\n', agcf, taufactor);
Daniel@0: 	[tar1 tar2 tar3 tar4; eps1 eps2 eps3 eps4]
Daniel@0: end
Daniel@0: 
Daniel@0: for i=0:nOutputSamples-1
Daniel@0: 	[sosOutput sosState]= soscascade(x(i*df+1:i*df+df), earFilters, ...  
Daniel@0: 				sosState);
Daniel@0: 	output = max(0, sosOutput);		%% Half Wave Rectify
Daniel@0: 	output(1) = 0;					%% Test Hack to make inversion easier.
Daniel@0: 	output(2) = 0;
Daniel@0: 	if agcf > 0
Daniel@0: 		[output agcState] = agc(output, [tar1 tar2 tar3 tar4; ...
Daniel@0: 				eps1 eps2 eps3 eps4], ...
Daniel@0: 				agcState);
Daniel@0: 	end
Daniel@0: 
Daniel@0: 	if differ > 0
Daniel@0: 		output = [output(1,:);output(1:nChannels-1,:) - ...
Daniel@0: 					output(2:nChannels,:)];
Daniel@0: 		output = max(0, output);
Daniel@0: 	end
Daniel@0: 
Daniel@0: 	if df > 1
Daniel@0: 		[output decState] = sosfilters(output, decFilt, decState);
Daniel@0: 	end
Daniel@0: 	y(:,i+1) = output(:,df);
Daniel@0: end
Daniel@0: 
Daniel@0: %y = min(y,2*tar4);
Daniel@0: y=y(3:nChannels,:);