Mercurial > hg > map

comparison testPrograms/testRP2.m @ 38:c2204b18f4a2 tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
End nov big change
author Ray Meddis <rmeddis@essex.ac.uk>
date Mon, 28 Nov 2011 13:34:28 +0000
parents
children
comparison
equal deleted inserted replaced
37:771a643d5c29 38:c2204b18f4a2
1 function testRP2(MAPparamsName,paramChanges)
2 % testIHC evaluates IHC I/O function
3 % multiple BF can be used but only one is easier to interpret.
4 % e.g. testRP(1000,'Normal',{});
5
6 global experiment method inputStimulusParams
7 global stimulusParameters IHC_VResp_VivoParams IHC_cilia_RPParams
8 savePath=path;
9 addpath (['..' filesep 'utilities'],['..' filesep 'MAP'])
10 dbstop if error
11
12 figure(4), clf,
13 set (gcf, 'name', ['IHC'])
14 set(gcf,'position',[613 354 360 322])
15 drawColors='rgbkmcy';
16 drawnow
17
18 if nargin<3
19 paramChanges=[];
20 end
21 if nargin<2
22 MAPparamsName='Normal';
23 end
24 if nargin<3
25 BF=800;
26 end
27
28 levels=-20:10:100;
29 nLevels=length(levels);
30 toneDuration=.05;
31 silenceDuration=.01;
32 sampleRate=50000;
33 dt=1/sampleRate;
34
35 allIHC_RP_peak=[];
36 allIHC_RP_dc=[];
37
38 %% Ruggero
39 %%Ruggero data
40 RuggeroData=[
41 0 2.00E-10;
42 10 5.00E-10;
43 20 1.50E-09;
44 30 2.50E-09;
45 40 5.30E-09;
46 50 1.00E-08;
47 60 1.70E-08;
48 70 2.50E-08;
49 80 4.00E-08;
50 90 6.00E-08;
51 100 1.50E-07;
52 110 3.00E-07;
53 ];
54
55
56 BF=10000;
57 targetFrequency=BF;
58
59 IHC_RP_peak=zeros(nLevels,1);
60 IHC_RP_min=zeros(nLevels,1);
61 IHC_RP_dc=zeros(nLevels,1);
62
63 time=dt:dt:toneDuration;
64
65 rampDuration=0.004;
66 rampTime=dt:dt:rampDuration;
67 ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
68 ones(1,length(time)-length(rampTime))];
69 ramp=ramp.*fliplr(ramp);
70
71 silence=zeros(1,round(silenceDuration/dt));
72
73 toneStartptr=length(silence)+1;
74 toneMidptr=toneStartptr+round(toneDuration/(2*dt)) -1;
75 toneEndptr=toneStartptr+round(toneDuration/dt) -1;
76
77 levelNo=0;
78 for leveldB=levels
79 levelNo=levelNo+1;
80 % replicate at all levels
81 amp=28e-6*10^(leveldB/20);
82
83 %% create signal (leveldB/ frequency)
84 inputSignal=amp*sin(2*pi*targetFrequency'*time);
85 inputSignal= ramp.*inputSignal;
86 inputSignal=[silence inputSignal silence];
87 inputStimulusParams.sampleRate=1/dt;
88 % global IHC_ciliaParams
89
90 %% disable efferent for fast processing
91
92 %% run the model
93 global DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV...
94 IHCoutput
95 AN_spikesOrProbability='probability';
96
97 MAP1_14(inputSignal, sampleRate, BF, ...
98 MAPparamsName, AN_spikesOrProbability, paramChanges);
99
100 % DRNL
101 DRNLoutput=DRNLoutput;
102 DRNL_peak(levelNo,1)=max(DRNLoutput(toneMidptr:toneEndptr));
103 DRNL_min(levelNo,1)=min(DRNLoutput(toneMidptr:toneEndptr));
104 DRNL_dc(levelNo,1)=mean(DRNLoutput(toneMidptr:toneEndptr));
105
106 end
107 %% plot DRNL
108 subplot(2,2,1)
109 referenceDisp=10e-9;
110 semilogy(levels,DRNL_peak, 'linewidth',2), hold on
111 semilogy(RuggeroData(:,1),RuggeroData(:,2),'o')
112 title(['BM: Ruggero ' num2str(BF) ' Hz'])
113 ylabel ('displacement(m)'), xlabel('dB SPL')
114 xlim([min(levels) max(levels)]), ylim([1e-10 1e-7])
115 grid on
116
117
118 %% Dallos
119 BF=800;
120 targetFrequency=BF;
121
122 IHC_RP_peak=zeros(nLevels,1);
123 IHC_RP_min=zeros(nLevels,1);
124 IHC_RP_dc=zeros(nLevels,1);
125
126 time=dt:dt:toneDuration;
127
128 rampDuration=0.004;
129 rampTime=dt:dt:rampDuration;
130 ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
131 ones(1,length(time)-length(rampTime))];
132 ramp=ramp.*fliplr(ramp);
133
134 silence=zeros(1,round(silenceDuration/dt));
135
136 toneStartptr=length(silence)+1;
137 toneMidptr=toneStartptr+round(toneDuration/(2*dt)) -1;
138 toneEndptr=toneStartptr+round(toneDuration/dt) -1;
139
140 levelNo=0;
141 for leveldB=levels
142 levelNo=levelNo+1;
143 % replicate at all levels
144 amp=28e-6*10^(leveldB/20);
145
146 %% create signal (leveldB/ frequency)
147 inputSignal=amp*sin(2*pi*targetFrequency'*time);
148 inputSignal= ramp.*inputSignal;
149 inputSignal=[silence inputSignal silence];
150 inputStimulusParams.sampleRate=1/dt;
151 % global IHC_ciliaParams
152
153 %% disable efferent for fast processing
154 method.DRNLSave=1;
155 method.IHC_cilia_RPSave=1;
156 method.IHCpreSynapseSave=1;
157 method.IHC_cilia_RPSave=1;
158 method.segmentDuration=-1;
159 moduleSequence=1:4;
160
161 %% run the model
162 global DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV...
163 IHCoutput
164 AN_spikesOrProbability='probability';
165
166 MAP1_14(inputSignal, sampleRate, BF, ...
167 MAPparamsName, AN_spikesOrProbability, paramChanges);
168
169 % DRNL
170 DRNLoutput=DRNLoutput;
171 DRNL_peak(levelNo,1)=max(DRNLoutput(toneMidptr:toneEndptr));
172 DRNL_min(levelNo,1)=min(DRNLoutput(toneMidptr:toneEndptr));
173 DRNL_dc(levelNo,1)=mean(DRNLoutput(toneMidptr:toneEndptr));
174
175 % cilia
176 IHC_ciliaData=IHC_cilia_output;
177 IHC_ciliaData=IHC_ciliaData;
178 IHC_cilia_peak(levelNo,1)=...
179 max(IHC_ciliaData(toneMidptr:toneEndptr));
180 IHC_cilia_min(levelNo,1)=...
181 min(IHC_ciliaData(toneMidptr:toneEndptr));
182 IHC_cilia_dc(levelNo,1)=...
183 mean(IHC_ciliaData(toneMidptr:toneEndptr));
184
185 % RP
186 IHC_RPData=IHCoutput;
187 IHC_RPData=IHC_RPData;
188 IHC_RP_peak(levelNo,1)=...
189 max(IHC_RPData(toneMidptr:toneEndptr));
190 IHC_RP_min(levelNo,1)=...
191 min(IHC_RPData(toneMidptr:toneEndptr));
192 IHC_RP_dc(levelNo,1)=...
193 mean(IHC_RPData(toneMidptr:toneEndptr));
194
195 end
196
197 %% Dallos and Harris data
198 %% plot receptor potentials
199 figure(4)
200 subplot(2,2,3)
201 % RP I/O function min and max
202 restingRP=IHC_RP_peak(1);
203 toPlot= [fliplr(IHC_RP_min(:,1)') IHC_RP_peak(:,1)'];
204 microPa= 28e-6*10.^(levels/20);
205 microPa=[-fliplr(microPa) microPa];
206 plot(microPa,toPlot, 'linewidth',2)
207
208 dallosx=[-0.9 -0.1 -0.001 0.001 0.01 0.9];
209 dallosy=[-8 -7.8 -6.5 11 16.5 22]/1000 + restingRP;
210 subplot(2,2,3)
211 hold on, plot(dallosx,dallosy, 'o')
212 plot([-1 1], [restingRP restingRP], 'r')
213 title(' Dallos(86) data at 800 Hz')
214 ylabel ('receptor potential(V)'), xlabel('Pa')
215 ylim([-0.08 -0.02]), xlim([-1 1])
216 grid on
217
218
219 %% Patuzzi
220 BF=7000;
221 targetFrequency=BF;
222
223 IHC_RP_peak=zeros(nLevels,1);
224 IHC_RP_min=zeros(nLevels,1);
225 IHC_RP_dc=zeros(nLevels,1);
226
227 time=dt:dt:toneDuration;
228
229 rampDuration=0.004;
230 rampTime=dt:dt:rampDuration;
231 ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
232 ones(1,length(time)-length(rampTime))];
233 ramp=ramp.*fliplr(ramp);
234
235 silence=zeros(1,round(silenceDuration/dt));
236
237 toneStartptr=length(silence)+1;
238 toneMidptr=toneStartptr+round(toneDuration/(2*dt)) -1;
239 toneEndptr=toneStartptr+round(toneDuration/dt) -1;
240
241 levelNo=0;
242 for leveldB=levels
243 levelNo=levelNo+1;
244 % replicate at all levels
245 amp=28e-6*10^(leveldB/20);
246
247 %% create signal (leveldB/ frequency)
248 inputSignal=amp*sin(2*pi*targetFrequency'*time);
249 inputSignal= ramp.*inputSignal;
250 inputSignal=[silence inputSignal silence];
251 inputStimulusParams.sampleRate=1/dt;
252 % global IHC_ciliaParams
253
254 %% run the model
255 global DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV...
256 IHCoutput
257 AN_spikesOrProbability='probability';
258
259 MAP1_14(inputSignal, sampleRate, BF, ...
260 MAPparamsName, AN_spikesOrProbability, paramChanges);
261
262 % RP
263 IHC_RPData=IHCoutput;
264 IHC_RP_peak(levelNo,1)=...
265 max(IHC_RPData(toneMidptr:toneEndptr));
266 IHC_RP_min(levelNo,1)=...
267 min(IHC_RPData(toneMidptr:toneEndptr));
268 IHC_RP_dc(levelNo,1)=...
269 mean(IHC_RPData(toneMidptr:toneEndptr));
270 end
271
272
273 %% RP I/O function min and max
274 figure(4)
275 subplot(2,2,4)
276 restingRP=IHC_RP_peak(1);
277 peakRP=max(IHC_RP_peak);
278 plot(levels, IHC_RP_peak, 'linewidth',2)
279 hold on
280 plot(levels, IHC_RP_dc, ':', 'linewidth',2)
281 plot([min(levels) max(levels)], [restingRP restingRP], 'r')
282 xlim([min(levels) max(levels)])
283 % animal data
284 sndLevel=[5 15 25 35 45 55 65 75];
285 RPanimal=restingRP+[0.5 2 4.6 5.8 6.4 7.2 8 10.2]/1000;
286 % could be misleading when restingRP changes
287 RPanimal=-0.060+[0.5 2 4.6 5.8 6.4 7.2 8 10.2]/1000;
288 hold on, plot(sndLevel,RPanimal,'o')
289
290 grid on
291 title(' 7 kHz Patuzzi')
292 ylabel ('RP(V) peak and DC'), xlabel('dB SPL')
293 ylim([-0.07 -0.04])
294 path(savePath);
295 disp(paramChanges)