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