-
Notifications
You must be signed in to change notification settings - Fork 7
/
model2018.m
274 lines (244 loc) · 10.1 KB
/
model2018.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
function output=model2018(sign,fs,fc,irregularities,storeflag,subject,sheraPo,IrrPct,non_linear_type,nH,nM,nL,clean,data_folder)
% function output=model2018(sign,fs,fc,irregularities,storeflag,subject,sheraPo,IrrPct,non_linear_type,nH,nM,nL,clean,data_folder)
%
% Computational model of the auditory periphery (Verhulst, Altoe, Vasilkov, 2018).
%
% Input:
% -sign: stimulus
% -fs: samplerate
% -fc: probe frequency or alternatively a string 'all' to probe all cochlear
% sections or 'half' to probe half sections, 'abr' to store the
% 401 sections used to compute the abr responses in Verhulst et al. 2017
% -irregularities: decide wether turn on (1) or off (0) irregularities and
% nonlinearities of the cochlear model (default 1 )
% -storeflag: string that sets what variables to store from the computation,
% each letter correspond to one desired output variable (e.g., 'avhl'
% to store acceleration, displacement, high and low spont. rate fibers.)
% See "output" for the corresponding characters. Default: 'avihmlebw'.
% -subject: number representing the seed to generate the random
% irregularities in the cochlear sections (default 1)
% -sheraPo: starting real part of the poles of the cochlear model
% it can be either an array with one value per BM section, or a
% single value far all section (default 0.6)
% -IrrPct: magnitude of random perturbations on the BM (irregularities, default 0.05=5%)
% -non_linear_type: select the type of nonlinearity in the BM model.
% Currently implemented
% 'vel'= instantaneouns nonlinearity based on local BM velocity (see Verhulst et al. 2012)
% 'none'= linear model
% -nH,nM,nL= number of high, medium and low spont. fibers employed to
% compute the response of cn and ic nuclei. Default 13,3,3. This
% parameters can be passed either as a single value for all sections or
% as an array with each value corresponding to a single CF location (and
% the length of those array must match the number of the section probed)
%
% Output:
% -output.v: BM velocity (store 'v')
% -output.y: BM displacement (store 'y')
% -output.emission: pressure output from the middle ear (store 'e')
% -output.cf: center frequencies (always stored)
% -output.ihc: IHC receptor potential (store 'i')
% -output.anfH: HSR fiber spike probability [0,1] (store 'h')
% -output.anfM: MSR fiber spike probability [0,1] (store 'm')
% -output.anfL: LSR fiber spike probability [0,1] (store 'l')
% -output.an_sum: summation of HSR, MSR and LSR per channel,i.e., the input to the cn model (storeflag 'b')
% -output.cn: cochlear nuclei output (storeflag 'b')
% -output.ic: IC (storeflag 'b')
% -output.w1,outupt.w3,output.w5: wave 1,3 and 5 (storeflag 'w')
% -output.fs_bm= sampling frequency of the bm simulations
% -output.fs_ihc= sample rate of the inner hair cell output
% -output.fs_an= sample rate of the an output
% -output.fs_abr= sample rate of the IC,CN and W1/3/5 outputs
% -clean: delete the temporary files generated by the python call (1) or
% not (0). Default 1
% -data_folder: string of folder where to save the temporary files
% generated by the python call (default current folder)
%
% References:
% Verhulst, S., Altoe, A., and Vasilkov, V. (2018). Computational
% modeling of the human auditory periphery: Auditorynerve responses,
% evoked potentials and hearing loss. Hearing Research, vol. 360, pp. 55-75.
%
% See also:
% Verhulst, S., Ernst, F., Garrett, M. and Vasilkov, V. (2018).
% Suprathreshold psychoacoustics and envelope-following response
% relations: Normal-hearing, synaptopathy and cochlear gain loss.
% Acta Acustica united with Acustica, vol. 104, pp. 800-803.
%
% -------------------------------------------------------------------------
% The model code and interface was written by Alessandro Altoe and Sarah
% Verhulst (copyright 2012,2014,2015,2016,2018) and is licensed under the
% UGent acadamic license (see details in license file that is part of this
% repository).
% The Verhulstetal2018Model consists of the following files:
% tridiag.so, cochlea_utils.c, run_model2018.py, model2018.m,
% cochlear_model2018.py, inner_hair_cell2018.py, auditory_nerve2017.py,
% ic_cn2017.py, ExampleSimulation.m, ExampleAnalysis.m, and the
% HI profiles in the Poles folder.
% -------------------------------------------------------------------------
DECIMATION=5;
[channels,idx]=min(size(sign));
if nargin<13
data_folder=[pwd(),'/'];
end
if nargin<13
clean=1;
end
if nargin < 12
nL=3;
end
if nargin < 11
nM = 3;
end
if nargin < 10
nH = 13;
end
if nargin < 9
non_linear_type = 'vel';
end
if nargin < 8
IrrPct = 0.05;
end
if nargin < 7
sheraPo = 0.06;
end
if nargin < 6
subject = 1;
end
if nargin < 5
storeflag = 'vihlmeb';
end
if nargin< 4
irregularities=ones(1,channels);
end
if clean == 1
fname2clean = []; %%% AO
end
if(numel(irregularities)==1)
irregularities=irregularities*ones(1,channels);
end
Fs=fs;
sectionsNo=1e3;
[channels,idx]=min(size(sign));
if(idx==2) %transpose it (python C-style row major order)
sign=sign';
end
stim=sign;
if(ischar(fc) && strcmp(fc,'all')) %if probing all sections 1001 output (1000 sections plus the middle ear)
l=sectionsNo;
elseif(ischar(fc) && strcmp(fc,'half')) %if probing half sections sections 1001 output (1000 sections plus the middle ear)
l=sectionsNo/2;
elseif(ischar(fc) && strcmp(fc,'abr')) %if probing half sections sections 1001 output (1000 sections plus the middle ear)
l=401;
else %else pass it as a column vector
[l,idx]=max(size(fc));
if(idx==2)
fc=fc';
end
fc=round(fc);
end
probes=fc;
act_path=pwd;
%%% Storing input.mat which will be read 'later' in Python
fname = 'input.mat';
save(fname,'stim','Fs','channels','subject','sheraPo','irregularities','probes',...
'sectionsNo','data_folder','storeflag','IrrPct','non_linear_type','nH','nM','nL','-v7');
%%% Running the model
[status,res] = system('python run_model2018.py','-echo');
switch status
case 0 %%% If the simulation succeeded then the stored results are read
cd(data_folder);
for i=1:channels
p=length(stim(i,:));
p2=ceil(length(stim(i,:))/DECIMATION);
fname = strcat('cf',int2str(i),'.mat');
tmp = load(fname,'cf');
output(i).cf = tmp.cf;
if clean; fname2clean{end+1} = fname; end
if strfind(storeflag,'v')
fname = strcat('v',int2str(i),'.mat');
tmp = load(fname,'Vsolution');
output(i).v = tmp.Vsolution;
if clean; fname2clean{end+1} = fname; end
end
if strfind(storeflag,'y')
fname = strcat('y',int2str(i),'.mat');
tmp = load(fname,'Ysolution');
output(i).y = tmp.Ysolution;
if clean; fname2clean{end+1} = fname; end
end
if strfind(storeflag,'e')
fname = strcat('emission',int2str(i),'.mat');
tmp = load(fname,'oto_emission');
output(i).e = tmp.oto_emission;
if clean; fname2clean{end+1} = fname; end
end
if strfind(storeflag,'i')
fname = strcat('ihc',int2str(i),'.mat');
tmp = load(fname,'Vm');
output(i).ihc = tmp.Vm;
if clean; fname2clean{end+1} = fname; end
end
if strfind(storeflag,'h')
fname = strcat('anfH',int2str(i),'.mat');
tmp = load(fname,'anfH');
output(i).anfH = tmp.anfH;
if clean; fname2clean{end+1} = fname; end
end
if strfind(storeflag,'m')
fname = strcat('anfM',int2str(i),'.mat');
tmp = load(fname,'anfM');
output(i).anfM = tmp.anfM;
if clean; fname2clean{end+1} = fname; end
end
if strfind(storeflag,'l')
fname = strcat('anfL',int2str(i),'.mat');
tmp = load(fname,'anfL');
output(i).anfL = tmp.anfL;
if clean; fname2clean{end+1} = fname; end
end
if strfind(storeflag,'b')
fname = strcat('cn',int2str(i),'.mat');
tmp = load(fname,'cn');
output(i).cn = tmp.cn;
if clean; fname2clean{end+1} = fname; end
fname = strcat('ic',int2str(i),'.mat');
tmp = load(fname,'ic');
output(i).ic = tmp.ic;
if clean; fname2clean{end+1} = fname; end
fname = strcat('AN',int2str(i),'.mat');
tmp = load(fname,'anSummed');
output(i).an_summed = tmp.anSummed;
if clean; fname2clean{end+1} = fname; end
end
if strfind(storeflag,'w')
fname = strcat('1w',int2str(i),'.mat');
tmp = load(fname,'w1');
output(i).w1 = tmp.w1;
if clean; fname2clean{end+1} = fname; end
fname = strcat('3w',int2str(i),'.mat');
tmp = load(fname,'w3');
output(i).w3 = tmp.w3;
if clean; fname2clean{end+1} = fname; end %%% AO
fname = strcat('5w',int2str(i),'.mat');
tmp = load(fname,'w5');
output(i).w5 = tmp.w5;
if clean; fname2clean{end+1} = fname; end %%% AO
end
output(i).fs_abr=fs/DECIMATION;
output(i).fs_an=fs/DECIMATION;
output(i).fs_bm=fs;
output(i).fs_ihc=fs;
end
if clean==1
for i = 1:length(fname2clean)
try
delete(fname2clean{i});
end
end
% delete *.mat
end
cd(act_path);
otherwise %%% Then the simulation did not succeed
disp(res) % it re-prints the Python warning...
error('Something went wrong (see the Python warning)')
end