view toolboxes/MIRtoolbox1.3.2/AuditoryToolbox/SeneffEar.m @ 0:e9a9cd732c1e tip

first hg version after svn

function y = SeneffEar(x, fs, plotChannel)

% Compute the response of Stephanie Seneff's Auditory Model.  The
% input waveform x (with sampling rate fs) is converted into fourty
% channels of auditory firing probabilities.
%
% This m-function is based on data from the following paper:
%	Benjamin D. Bryant and John D. Gowdy, "Simulation of Stages
%	I and II of Seneff's Auditory Model (SAM) Using Matlab", and 
%	published in the Proceedings of the 1993 Matlab User's Group
%	Conference.
% Thanks to Benjamin Bryant for supplying us with his filter 
% coefficients and the initial organization of this implementation.

% (c) 1998 Interval Research Corporation  

global SeneffFS SeneffPreemphasis SeneffFilterBank
global SeneffForward SeneffBackward

if nargin < 2; fs = 16000; end
if nargin < 3; plotChannel = 0; end

if exist('SeneffFS') ~= 1 | length(SeneffFS) == 0 | SeneffFS ~= fs
	[SeneffPreemphasis, SeneffFilterBank, ...
	 SeneffForward, SeneffBackward] = SeneffEarSetup(fs);
	SeneffFS = fs;
end

y=soscascade(filter(SeneffPreemphasis, [1], x), ...
			SeneffFilterBank);

[width,channels] = size(SeneffForward);
for j=1:channels
	y(j,:) = filter(SeneffForward(:,j),SeneffBackward(:,j),y(j,:));
end

if plotChannel > 0
	clf
	subplot(4,2,1);
	plot(y(plotChannel,:))
	subplot(4,2,2);
	f=find(y(plotChannel,:)>max(y(plotChannel,:))/20);
	plotStart = max(1,f(1)-10);
	plotEnd = min(plotStart+84,length(y(plotChannel,:)));
	plot(y(plotChannel,plotStart:plotEnd));
	ylabel('Filter Bank')
end

% Implement Seneff's detector non-linearity.
hwrA = 10;
hwrB = 65;
hwrG = 2.35;

y = hwrA*atan(hwrB*max(0,y))+exp(hwrA*hwrB*min(0,y));
if plotChannel > 0
	subplot(4,2,3);
	plot(y(plotChannel,:))
	subplot(4,2,4);
	plot(y(plotChannel,plotStart:plotEnd));
	ylabel('HWR');
end

% Implement Seneff's short-term adaptation (a reservoir hair cell
% model.)
start_Cn = 442.96875;
Tua = 58.3333/fs;
Tub = 8.3333/fs;

len = length(x);
Cn = start_Cn*ones(channels,1)*0;		% Does this work non zero?
for j=1:len
	Sn = y(:,j);
	flow = max(0,Tua*(Sn-Cn));
	Cn = Cn + flow - Tub*Cn;
	y(:,j) = flow;
end

% Implement Seneff's Low Pass filter (to model the loss of
% Synchrony.

lpAlpha = .209611;
lpPoly = poly([lpAlpha lpAlpha lpAlpha lpAlpha]);
Glp = sum(lpPoly);

for j=1:channels
	y(j,:) = filter([Glp],lpPoly,y(j,:));
end

if plotChannel > 0
	subplot(4,2,5);
	plot(y(plotChannel,:))
	subplot(4,2,6);
	plot(y(plotChannel,plotStart:plotEnd));
	ylabel('Adaptation');
end

% Finally implement Seneff's Adaptive Gain Control.  This computes
%				y(n) = y(n)/(1+k y'(n))
% where y'(n) is a low-pass filtered version of the input.  Note, 
% this is a single channel AGC.

initial_yn = 0.23071276;
alpha_agc = 0.979382181;
kagc = 0.002;

for j=1:channels
	averageY(j,:) = filter([0 1-alpha_agc],[alpha_agc],y(j,:),initial_yn);
end
y = y./(1+kagc*averageY);

if plotChannel > 0
	subplot(4,2,7);
	plot(y(plotChannel,:))
	subplot(4,2,8);
	plot(y(plotChannel,plotStart:plotEnd));
	ylabel('AGC');
end