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