wolffd@0: function y = SeneffEar(x, fs, plotChannel)
wolffd@0: 
wolffd@0: % Compute the response of Stephanie Seneff's Auditory Model.  The
wolffd@0: % input waveform x (with sampling rate fs) is converted into fourty
wolffd@0: % channels of auditory firing probabilities.
wolffd@0: %
wolffd@0: % This m-function is based on data from the following paper:
wolffd@0: %	Benjamin D. Bryant and John D. Gowdy, "Simulation of Stages
wolffd@0: %	I and II of Seneff's Auditory Model (SAM) Using Matlab", and 
wolffd@0: %	published in the Proceedings of the 1993 Matlab User's Group
wolffd@0: %	Conference.
wolffd@0: % Thanks to Benjamin Bryant for supplying us with his filter 
wolffd@0: % coefficients and the initial organization of this implementation.
wolffd@0: 
wolffd@0: % (c) 1998 Interval Research Corporation  
wolffd@0: 
wolffd@0: global SeneffFS SeneffPreemphasis SeneffFilterBank
wolffd@0: global SeneffForward SeneffBackward
wolffd@0: 
wolffd@0: if nargin < 2; fs = 16000; end
wolffd@0: if nargin < 3; plotChannel = 0; end
wolffd@0: 
wolffd@0: if exist('SeneffFS') ~= 1 | length(SeneffFS) == 0 | SeneffFS ~= fs
wolffd@0: 	[SeneffPreemphasis, SeneffFilterBank, ...
wolffd@0: 	 SeneffForward, SeneffBackward] = SeneffEarSetup(fs);
wolffd@0: 	SeneffFS = fs;
wolffd@0: end
wolffd@0: 
wolffd@0: y=soscascade(filter(SeneffPreemphasis, [1], x), ...
wolffd@0: 			SeneffFilterBank);
wolffd@0: 
wolffd@0: [width,channels] = size(SeneffForward);
wolffd@0: for j=1:channels
wolffd@0: 	y(j,:) = filter(SeneffForward(:,j),SeneffBackward(:,j),y(j,:));
wolffd@0: end
wolffd@0: 
wolffd@0: if plotChannel > 0
wolffd@0: 	clf
wolffd@0: 	subplot(4,2,1);
wolffd@0: 	plot(y(plotChannel,:))
wolffd@0: 	subplot(4,2,2);
wolffd@0: 	f=find(y(plotChannel,:)>max(y(plotChannel,:))/20);
wolffd@0: 	plotStart = max(1,f(1)-10);
wolffd@0: 	plotEnd = min(plotStart+84,length(y(plotChannel,:)));
wolffd@0: 	plot(y(plotChannel,plotStart:plotEnd));
wolffd@0: 	ylabel('Filter Bank')
wolffd@0: end
wolffd@0: 
wolffd@0: % Implement Seneff's detector non-linearity.
wolffd@0: hwrA = 10;
wolffd@0: hwrB = 65;
wolffd@0: hwrG = 2.35;
wolffd@0: 
wolffd@0: y = hwrA*atan(hwrB*max(0,y))+exp(hwrA*hwrB*min(0,y));
wolffd@0: if plotChannel > 0
wolffd@0: 	subplot(4,2,3);
wolffd@0: 	plot(y(plotChannel,:))
wolffd@0: 	subplot(4,2,4);
wolffd@0: 	plot(y(plotChannel,plotStart:plotEnd));
wolffd@0: 	ylabel('HWR');
wolffd@0: end
wolffd@0: 
wolffd@0: % Implement Seneff's short-term adaptation (a reservoir hair cell
wolffd@0: % model.)
wolffd@0: start_Cn = 442.96875;
wolffd@0: Tua = 58.3333/fs;
wolffd@0: Tub = 8.3333/fs;
wolffd@0: 
wolffd@0: len = length(x);
wolffd@0: Cn = start_Cn*ones(channels,1)*0;		% Does this work non zero?
wolffd@0: for j=1:len
wolffd@0: 	Sn = y(:,j);
wolffd@0: 	flow = max(0,Tua*(Sn-Cn));
wolffd@0: 	Cn = Cn + flow - Tub*Cn;
wolffd@0: 	y(:,j) = flow;
wolffd@0: end
wolffd@0: 
wolffd@0: % Implement Seneff's Low Pass filter (to model the loss of
wolffd@0: % Synchrony.
wolffd@0: 
wolffd@0: lpAlpha = .209611;
wolffd@0: lpPoly = poly([lpAlpha lpAlpha lpAlpha lpAlpha]);
wolffd@0: Glp = sum(lpPoly);
wolffd@0: 
wolffd@0: for j=1:channels
wolffd@0: 	y(j,:) = filter([Glp],lpPoly,y(j,:));
wolffd@0: end
wolffd@0: 
wolffd@0: if plotChannel > 0
wolffd@0: 	subplot(4,2,5);
wolffd@0: 	plot(y(plotChannel,:))
wolffd@0: 	subplot(4,2,6);
wolffd@0: 	plot(y(plotChannel,plotStart:plotEnd));
wolffd@0: 	ylabel('Adaptation');
wolffd@0: end
wolffd@0: 
wolffd@0: % Finally implement Seneff's Adaptive Gain Control.  This computes
wolffd@0: %				y(n) = y(n)/(1+k y'(n))
wolffd@0: % where y'(n) is a low-pass filtered version of the input.  Note, 
wolffd@0: % this is a single channel AGC.
wolffd@0: 
wolffd@0: initial_yn = 0.23071276;
wolffd@0: alpha_agc = 0.979382181;
wolffd@0: kagc = 0.002;
wolffd@0: 
wolffd@0: for j=1:channels
wolffd@0: 	averageY(j,:) = filter([0 1-alpha_agc],[alpha_agc],y(j,:),initial_yn);
wolffd@0: end
wolffd@0: y = y./(1+kagc*averageY);
wolffd@0: 
wolffd@0: if plotChannel > 0
wolffd@0: 	subplot(4,2,7);
wolffd@0: 	plot(y(plotChannel,:))
wolffd@0: 	subplot(4,2,8);
wolffd@0: 	plot(y(plotChannel,plotStart:plotEnd));
wolffd@0: 	ylabel('AGC');
wolffd@0: end
wolffd@0: