|
wolffd@0
|
1 function y = SeneffEar(x, fs, plotChannel)
|
|
wolffd@0
|
2
|
|
wolffd@0
|
3 % Compute the response of Stephanie Seneff's Auditory Model. The
|
|
wolffd@0
|
4 % input waveform x (with sampling rate fs) is converted into fourty
|
|
wolffd@0
|
5 % channels of auditory firing probabilities.
|
|
wolffd@0
|
6 %
|
|
wolffd@0
|
7 % This m-function is based on data from the following paper:
|
|
wolffd@0
|
8 % Benjamin D. Bryant and John D. Gowdy, "Simulation of Stages
|
|
wolffd@0
|
9 % I and II of Seneff's Auditory Model (SAM) Using Matlab", and
|
|
wolffd@0
|
10 % published in the Proceedings of the 1993 Matlab User's Group
|
|
wolffd@0
|
11 % Conference.
|
|
wolffd@0
|
12 % Thanks to Benjamin Bryant for supplying us with his filter
|
|
wolffd@0
|
13 % coefficients and the initial organization of this implementation.
|
|
wolffd@0
|
14
|
|
wolffd@0
|
15 % (c) 1998 Interval Research Corporation
|
|
wolffd@0
|
16
|
|
wolffd@0
|
17 global SeneffFS SeneffPreemphasis SeneffFilterBank
|
|
wolffd@0
|
18 global SeneffForward SeneffBackward
|
|
wolffd@0
|
19
|
|
wolffd@0
|
20 if nargin < 2; fs = 16000; end
|
|
wolffd@0
|
21 if nargin < 3; plotChannel = 0; end
|
|
wolffd@0
|
22
|
|
wolffd@0
|
23 if exist('SeneffFS') ~= 1 | length(SeneffFS) == 0 | SeneffFS ~= fs
|
|
wolffd@0
|
24 [SeneffPreemphasis, SeneffFilterBank, ...
|
|
wolffd@0
|
25 SeneffForward, SeneffBackward] = SeneffEarSetup(fs);
|
|
wolffd@0
|
26 SeneffFS = fs;
|
|
wolffd@0
|
27 end
|
|
wolffd@0
|
28
|
|
wolffd@0
|
29 y=soscascade(filter(SeneffPreemphasis, [1], x), ...
|
|
wolffd@0
|
30 SeneffFilterBank);
|
|
wolffd@0
|
31
|
|
wolffd@0
|
32 [width,channels] = size(SeneffForward);
|
|
wolffd@0
|
33 for j=1:channels
|
|
wolffd@0
|
34 y(j,:) = filter(SeneffForward(:,j),SeneffBackward(:,j),y(j,:));
|
|
wolffd@0
|
35 end
|
|
wolffd@0
|
36
|
|
wolffd@0
|
37 if plotChannel > 0
|
|
wolffd@0
|
38 clf
|
|
wolffd@0
|
39 subplot(4,2,1);
|
|
wolffd@0
|
40 plot(y(plotChannel,:))
|
|
wolffd@0
|
41 subplot(4,2,2);
|
|
wolffd@0
|
42 f=find(y(plotChannel,:)>max(y(plotChannel,:))/20);
|
|
wolffd@0
|
43 plotStart = max(1,f(1)-10);
|
|
wolffd@0
|
44 plotEnd = min(plotStart+84,length(y(plotChannel,:)));
|
|
wolffd@0
|
45 plot(y(plotChannel,plotStart:plotEnd));
|
|
wolffd@0
|
46 ylabel('Filter Bank')
|
|
wolffd@0
|
47 end
|
|
wolffd@0
|
48
|
|
wolffd@0
|
49 % Implement Seneff's detector non-linearity.
|
|
wolffd@0
|
50 hwrA = 10;
|
|
wolffd@0
|
51 hwrB = 65;
|
|
wolffd@0
|
52 hwrG = 2.35;
|
|
wolffd@0
|
53
|
|
wolffd@0
|
54 y = hwrA*atan(hwrB*max(0,y))+exp(hwrA*hwrB*min(0,y));
|
|
wolffd@0
|
55 if plotChannel > 0
|
|
wolffd@0
|
56 subplot(4,2,3);
|
|
wolffd@0
|
57 plot(y(plotChannel,:))
|
|
wolffd@0
|
58 subplot(4,2,4);
|
|
wolffd@0
|
59 plot(y(plotChannel,plotStart:plotEnd));
|
|
wolffd@0
|
60 ylabel('HWR');
|
|
wolffd@0
|
61 end
|
|
wolffd@0
|
62
|
|
wolffd@0
|
63 % Implement Seneff's short-term adaptation (a reservoir hair cell
|
|
wolffd@0
|
64 % model.)
|
|
wolffd@0
|
65 start_Cn = 442.96875;
|
|
wolffd@0
|
66 Tua = 58.3333/fs;
|
|
wolffd@0
|
67 Tub = 8.3333/fs;
|
|
wolffd@0
|
68
|
|
wolffd@0
|
69 len = length(x);
|
|
wolffd@0
|
70 Cn = start_Cn*ones(channels,1)*0; % Does this work non zero?
|
|
wolffd@0
|
71 for j=1:len
|
|
wolffd@0
|
72 Sn = y(:,j);
|
|
wolffd@0
|
73 flow = max(0,Tua*(Sn-Cn));
|
|
wolffd@0
|
74 Cn = Cn + flow - Tub*Cn;
|
|
wolffd@0
|
75 y(:,j) = flow;
|
|
wolffd@0
|
76 end
|
|
wolffd@0
|
77
|
|
wolffd@0
|
78 % Implement Seneff's Low Pass filter (to model the loss of
|
|
wolffd@0
|
79 % Synchrony.
|
|
wolffd@0
|
80
|
|
wolffd@0
|
81 lpAlpha = .209611;
|
|
wolffd@0
|
82 lpPoly = poly([lpAlpha lpAlpha lpAlpha lpAlpha]);
|
|
wolffd@0
|
83 Glp = sum(lpPoly);
|
|
wolffd@0
|
84
|
|
wolffd@0
|
85 for j=1:channels
|
|
wolffd@0
|
86 y(j,:) = filter([Glp],lpPoly,y(j,:));
|
|
wolffd@0
|
87 end
|
|
wolffd@0
|
88
|
|
wolffd@0
|
89 if plotChannel > 0
|
|
wolffd@0
|
90 subplot(4,2,5);
|
|
wolffd@0
|
91 plot(y(plotChannel,:))
|
|
wolffd@0
|
92 subplot(4,2,6);
|
|
wolffd@0
|
93 plot(y(plotChannel,plotStart:plotEnd));
|
|
wolffd@0
|
94 ylabel('Adaptation');
|
|
wolffd@0
|
95 end
|
|
wolffd@0
|
96
|
|
wolffd@0
|
97 % Finally implement Seneff's Adaptive Gain Control. This computes
|
|
wolffd@0
|
98 % y(n) = y(n)/(1+k y'(n))
|
|
wolffd@0
|
99 % where y'(n) is a low-pass filtered version of the input. Note,
|
|
wolffd@0
|
100 % this is a single channel AGC.
|
|
wolffd@0
|
101
|
|
wolffd@0
|
102 initial_yn = 0.23071276;
|
|
wolffd@0
|
103 alpha_agc = 0.979382181;
|
|
wolffd@0
|
104 kagc = 0.002;
|
|
wolffd@0
|
105
|
|
wolffd@0
|
106 for j=1:channels
|
|
wolffd@0
|
107 averageY(j,:) = filter([0 1-alpha_agc],[alpha_agc],y(j,:),initial_yn);
|
|
wolffd@0
|
108 end
|
|
wolffd@0
|
109 y = y./(1+kagc*averageY);
|
|
wolffd@0
|
110
|
|
wolffd@0
|
111 if plotChannel > 0
|
|
wolffd@0
|
112 subplot(4,2,7);
|
|
wolffd@0
|
113 plot(y(plotChannel,:))
|
|
wolffd@0
|
114 subplot(4,2,8);
|
|
wolffd@0
|
115 plot(y(plotChannel,plotStart:plotEnd));
|
|
wolffd@0
|
116 ylabel('AGC');
|
|
wolffd@0
|
117 end
|
|
wolffd@0
|
118
|
