Mercurial > hg > map
changeset 9:ecad0ea62b43
May27 mainly better parameters
| author | Ray Meddis <rmeddis@essex.ac.uk> |
|---|---|
| date | Tue, 31 May 2011 09:13:07 +0100 |
| parents | eafe11c86f44 |
| children | 4ef8e0e63866 |
| files | MAP/MAP1_14.m MAP/MAP1_14parallel.m multithreshold 1.46/expGUI_MT.fig multithreshold 1.46/expGUI_MT.m multithreshold 1.46/hs_err_pid6852.log multithreshold 1.46/old files/MAPmodel.m multithreshold 1.46/paradigms/paradigm_thr_IFMC.m multithreshold 1.46/paradigms/paradigm_thr_TMC.m multithreshold 1.46/paradigms/paradigm_training.m multithreshold 1.46/printReport.m multithreshold 1.46/savedData/mostRecentResults.mat multithreshold 1.46/subjGUI_MT.m multithreshold 1.46/testAN.m multithreshold 1.46/testBM.m multithreshold 1.46/testPhaseLocking.m multithreshold 1.46/testRP.m multithreshold 1.46/testSynapse.m parameterStore/MAPparamsNormal.m testPrograms/demoTwisterProbability.m testPrograms/demoTwisterSpikes.m testPrograms/test_MAP1_14.m |
| diffstat | 21 files changed, 2319 insertions(+), 141 deletions(-) [+] |
line wrap: on
line diff
--- a/MAP/MAP1_14.m Fri May 27 17:50:40 2011 +0100 +++ b/MAP/MAP1_14.m Tue May 31 09:13:07 2011 +0100 @@ -298,8 +298,7 @@ IHCGa= IHC_cilia_RPParams.Ga; % (leakage) IHCGu0 = IHCGa+IHCGmax./(1+exp(IHCu0/IHCs0).*(1+exp(IHCu1/IHCs1))); - - +IHCrestingCiliaCond=IHCGu0; % Receptor potential IHC_Cab= IHC_cilia_RPParams.Cab; @@ -308,9 +307,7 @@ IHC_Ek= IHC_cilia_RPParams.Ek; IHC_Ekp= IHC_Ek+IHC_Et*IHC_cilia_RPParams.Rpc; -IHCrestingV= -0.06; - -IHCrestingV= (IHC_Gk*IHC_Ek+IHCGu0*IHC_Et)/(IHCGu0+IHC_Gk); +IHCrestingV= (IHC_Gk*IHC_Ekp+IHCGu0*IHC_Et)/(IHCGu0+IHC_Gk); IHC_Vnow= IHCrestingV*ones(nBFs,1); % initial voltage IHC_RP= zeros(nBFs,segmentLength); @@ -677,9 +674,9 @@ IHCciliaDisplacement; % compute apical conductance - G=1./(1+exp(-(IHCciliaDisplacement-IHCu0)/IHCs0).*... + G=IHCGmax./(1+exp(-(IHCciliaDisplacement-IHCu0)/IHCs0).*... (1+exp(-(IHCciliaDisplacement-IHCu1)/IHCs1))); - Gu=IHCGmax*G + IHCGa; + Gu=G + IHCGa; % Compute receptor potential for idx=1:segmentLength @@ -1106,4 +1103,4 @@ end % segment -path(restorePath) \ No newline at end of file +path(restorePath)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/MAP/MAP1_14parallel.m Tue May 31 09:13:07 2011 +0100 @@ -0,0 +1,1148 @@ + +function MAP1_14parallel(inputSignal, sampleRate, BFlist, MAPparamsName, ... + AN_spikesOrProbability, paramChanges) +% To test this function use test_MAP1_14 in this folder +% +% All arguments are mandatory. +% +% BFlist is a list of BFs but can be '-1' to allow MAPparams to choose +% + +% MAPparamsName='Normal'; % source of model parameters +% AN_spikesOrProbability='spikes'; % or 'probability' +% paramChanges is a cell array of strings that can be used to make last +% minute parameter changes, e.g., to simulate OHC loss +% paramChanges{1}= 'DRNLParams.a=0;'; + +% The model parameters are established in the MAPparams<***> file +% and stored as global + +restorePath=path; +addpath (['..' filesep 'parameterStore']) + +global OMEParams DRNLParams IHC_cilia_RPParams IHCpreSynapseParams +global AN_IHCsynapseParams MacGregorParams MacGregorMultiParams + +% All of the results of this function are stored as global +global dt ANdt savedBFlist saveAN_spikesOrProbability saveMAPparamsName... + savedInputSignal OMEextEarPressure TMoutput OMEoutput ARattenuation ... + DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV... + IHCoutput ANprobRateOutput ANoutput savePavailable tauCas ... + CNoutput ICoutput ICmembraneOutput ICfiberTypeRates MOCattenuation + +% Normally only ICoutput(logical spike matrix) or ANprobRateOutput will be +% needed by the user; so the following will suffice +% global ANdt ICoutput ANprobRateOutput + +% Note that sampleRate has not changed from the original function call and +% ANprobRateOutput is sampled at this rate +% However ANoutput, CNoutput and IC output are stored as logical +% 'spike' matrices using a lower sample rate (see ANdt). + +% When AN_spikesOrProbability is set to probability, +% no spike matrices are computed. +% When AN_spikesOrProbability is set to 'spikes', +% no probability output is computed + +% Efferent control variables are ARattenuation and MOCattenuation +% These are scalars between 1 (no attenuation) and 0. +% They are represented with dt=1/sampleRate (not ANdt) +% They are computed using either AN probability rate output +% or IC (spikes) output as approrpriate. +% AR is computed using across channel activity +% MOC is computed on a within-channel basis. + + +% save as global for later plotting if required +savedBFlist=BFlist; +saveAN_spikesOrProbability=AN_spikesOrProbability; +saveMAPparamsName=MAPparamsName; + +% Read parameters from MAPparams<***> file in 'parameterStore' folder +cmd=['method=MAPparams' MAPparamsName ... + '(BFlist, sampleRate, 0);']; +eval(cmd); + +% Beware, 'BFlist=-1' is a legitimate argument for MAPparams<> +% if the calling program allows MAPparams to specify the list +BFlist=DRNLParams.nonlinCFs; + +% now accept last mintue parameter changes required by the calling program +if nargin>5 && ~isempty(paramChanges) + nChanges=length(paramChanges); + for idx=1:nChanges + eval(paramChanges{idx}) + end +end + +dt=1/sampleRate; +duration=length(inputSignal)/sampleRate; +% segmentDuration is specified in parameter file (must be >efferent delay) +segmentDuration=method.segmentDuration; +segmentLength=round(segmentDuration/ dt); +segmentTime=dt*(1:segmentLength); % used in debugging plots + +% all spiking activity is computed using longer epochs +ANspeedUpFactor=5; % 5 times longer + +% inputSignal must be row vector +[r c]=size(inputSignal); +if r>c, inputSignal=inputSignal'; end % transpose +% ignore stereo signals +inputSignal=inputSignal(1,:); % drop any second channel +savedInputSignal=inputSignal; + +% Segment the signal +% The sgment length is given but the signal length must be adjusted to be a +% multiple of both the segment length and the reduced segmentlength +[nSignalRows signalLength]=size(inputSignal); +segmentLength=ceil(segmentLength/ANspeedUpFactor)*ANspeedUpFactor; +% Make the signal length a whole multiple of the segment length +nSignalSegments=ceil(signalLength/segmentLength); +padSize=nSignalSegments*segmentLength-signalLength; +pad=zeros(nSignalRows,padSize); +inputSignal=[inputSignal pad]; +[ignore signalLength]=size(inputSignal); + +% AN (spikes) is computed at a lower sample rate when spikes required +% so it has a reduced segment length (see 'ANspeeUpFactor' above) +% AN CN and IC all use this sample interval +ANdt=dt*ANspeedUpFactor; +reducedSegmentLength=round(segmentLength/ANspeedUpFactor); +reducedSignalLength= round(signalLength/ANspeedUpFactor); + +%% Initialise with respect to each stage before computing +% by allocating memory, +% by computing constants +% by establishing easy to read variable names +% The computations are made in segments and boundary conditions must +% be established and stored. These are found in variables with +% 'boundary' or 'bndry' in the name + +%% OME --- +% external ear resonances +OMEexternalResonanceFilters=OMEParams.externalResonanceFilters; +[nOMEExtFilters c]=size(OMEexternalResonanceFilters); +% details of external (outer ear) resonances +OMEgaindBs=OMEexternalResonanceFilters(:,1); +OMEgainScalars=10.^(OMEgaindBs/20); +OMEfilterOrder=OMEexternalResonanceFilters(:,2); +OMElowerCutOff=OMEexternalResonanceFilters(:,3); +OMEupperCutOff=OMEexternalResonanceFilters(:,4); +% external resonance coefficients +ExtFilter_b=cell(nOMEExtFilters,1); +ExtFilter_a=cell(nOMEExtFilters,1); +for idx=1:nOMEExtFilters + Nyquist=sampleRate/2; + [b, a] = butter(OMEfilterOrder(idx), ... + [OMElowerCutOff(idx) OMEupperCutOff(idx)]... + /Nyquist); + ExtFilter_b{idx}=b; + ExtFilter_a{idx}=a; +end +OMEExtFilterBndry=cell(2,1); +OMEextEarPressure=zeros(1,signalLength); % pressure at tympanic membrane + +% pressure to velocity conversion using smoothing filter (50 Hz cutoff) +tau=1/(2*pi*50); +a1=dt/tau-1; a0=1; +b0=1+ a1; +TMdisp_b=b0; TMdisp_a=[a0 a1]; +% figure(9), freqz(TMdisp_b, TMdisp_a) +OME_TMdisplacementBndry=[]; + +% OME high pass (simulates poor low frequency stapes response) +OMEhighPassHighCutOff=OMEParams.OMEstapesLPcutoff; +Nyquist=sampleRate/2; +[stapesDisp_b,stapesDisp_a] = butter(1, OMEhighPassHighCutOff/Nyquist, 'high'); +% figure(10), freqz(stapesDisp_b, stapesDisp_a) + +OMEhighPassBndry=[]; + +% OMEampStapes might be reducdant (use OMEParams.stapesScalar) +stapesScalar= OMEParams.stapesScalar; + +% Acoustic reflex +efferentDelayPts=round(OMEParams.ARdelay/dt); +% smoothing filter +% Nyquist=(1/ANdt)/2; +% [ARfilt_b,ARfilt_a] = butter(1, (1/(2*pi*OMEParams.ARtau))/Nyquist, 'low'); +a1=dt/OMEParams.ARtau-1; a0=1; +b0=1+ a1; +ARfilt_b=b0; ARfilt_a=[a0 a1]; + +ARattenuation=ones(1,signalLength); +ARrateThreshold=OMEParams.ARrateThreshold; % may not be used +ARrateToAttenuationFactor=OMEParams.rateToAttenuationFactor; +ARrateToAttenuationFactorProb=OMEParams.rateToAttenuationFactorProb; +ARboundary=[]; +ARboundaryProb=0; + +% save complete OME record (stapes displacement) +OMEoutput=zeros(1,signalLength); +TMoutput=zeros(1,signalLength); + +%% BM --- +% BM is represented as a list of locations identified by BF +DRNL_BFs=BFlist; +nBFs= length(DRNL_BFs); + +% DRNLchannelParameters=DRNLParams.channelParameters; +DRNLresponse= zeros(nBFs, segmentLength); + +MOCrateToAttenuationFactor=DRNLParams.rateToAttenuationFactor; +rateToAttenuationFactorProb=DRNLParams.rateToAttenuationFactorProb; +MOCrateThreshold=DRNLParams.MOCrateThreshold; + +% smoothing filter for MOC +% Nyquist=(1/ANdt)/2; +% [MOCfilt_b,MOCfilt_a] = ... +% butter(1, (1/(2*pi*DRNLParams.MOCtau))/Nyquist, 'low'); +% figure(10), freqz(stapesDisp_b, stapesDisp_a) +a1=dt/DRNLParams.MOCtau-1; a0=1; +b0=1+ a1; +MOCfilt_b=b0; MOCfilt_a=[a0 a1]; +% figure(9), freqz(stapesDisp_b, stapesDisp_a) +MOCboundary=cell(nBFs,1); +MOCprobBoundary=cell(nBFs,1); + +MOCattSegment=zeros(nBFs,reducedSegmentLength); +MOCattenuation=ones(nBFs,signalLength); + +if DRNLParams.a>0 + DRNLcompressionThreshold=10^((1/(1-DRNLParams.c))* ... + log10(DRNLParams.b/DRNLParams.a)); +else + DRNLcompressionThreshold=inf; +end + +DRNLlinearOrder= DRNLParams.linOrder; +DRNLnonlinearOrder= DRNLParams.nonlinOrder; + +DRNLa=DRNLParams.a; +DRNLb=DRNLParams.b; +DRNLc=DRNLParams.c; +linGAIN=DRNLParams.g; +% +% gammatone filter coefficients for linear pathway +bw=DRNLParams.linBWs'; +phi = 2 * pi * bw * dt; +cf=DRNLParams.linCFs'; +theta = 2 * pi * cf * dt; +cos_theta = cos(theta); +sin_theta = sin(theta); +alpha = -exp(-phi).* cos_theta; +b0 = ones(nBFs,1); +b1 = 2 * alpha; +b2 = exp(-2 * phi); +z1 = (1 + alpha .* cos_theta) - (alpha .* sin_theta) * i; +z2 = (1 + b1 .* cos_theta) - (b1 .* sin_theta) * i; +z3 = (b2 .* cos(2 * theta)) - (b2 .* sin(2 * theta)) * i; +tf = (z2 + z3) ./ z1; +a0 = abs(tf); +a1 = alpha .* a0; +GTlin_a = [b0, b1, b2]; +GTlin_b = [a0, a1]; +GTlinOrder=DRNLlinearOrder; +GTlinBdry1=cell(nBFs); +GTlinBdry2=cell(nBFs); +GTlinBdry3=cell(nBFs); +GTlinBdry1out=cell(nBFs); +GTlinBdry2out=cell(nBFs); +GTlinBdry3out=cell(nBFs); + +% nonlinear gammatone filter coefficients +bw=DRNLParams.nlBWs'; +phi = 2 * pi * bw * dt; +cf=DRNLParams.nonlinCFs'; +theta = 2 * pi * cf * dt; +cos_theta = cos(theta); +sin_theta = sin(theta); +alpha = -exp(-phi).* cos_theta; +b0 = ones(nBFs,1); +b1 = 2 * alpha; +b2 = exp(-2 * phi); +z1 = (1 + alpha .* cos_theta) - (alpha .* sin_theta) * i; +z2 = (1 + b1 .* cos_theta) - (b1 .* sin_theta) * i; +z3 = (b2 .* cos(2 * theta)) - (b2 .* sin(2 * theta)) * i; +tf = (z2 + z3) ./ z1; +a0 = abs(tf); +a1 = alpha .* a0; +GTnonlin_a = [b0, b1, b2]; +GTnonlin_b = [a0, a1]; +GTnonlinOrder=DRNLnonlinearOrder; +firstGTnonlinBdry1=cell(nBFs); +firstGTnonlinBdry2=cell(nBFs); +firstGTnonlinBdry3=cell(nBFs); +firstGTnonlinBdry1out=cell(nBFs); +firstGTnonlinBdry2out=cell(nBFs); +firstGTnonlinBdry3out=cell(nBFs); + +secondGTnonlinBdry1=cell(nBFs); +secondGTnonlinBdry2=cell(nBFs); +secondGTnonlinBdry3=cell(nBFs); +secondGTnonlinBdry1out=cell(nBFs); +secondGTnonlinBdry2out=cell(nBFs); +secondGTnonlinBdry3out=cell(nBFs); + +% complete BM record (BM displacement) +DRNLoutput=zeros(nBFs, signalLength); + + +%% IHC --- +% IHC cilia activity and receptor potential +% viscous coupling between BM and stereocilia displacement +% Nyquist=sampleRate/2; +% IHCcutoff=1/(2*pi*IHC_cilia_RPParams.tc); +% [IHCciliaFilter_b,IHCciliaFilter_a]=... +% butter(1, IHCcutoff/Nyquist, 'high'); +a1=dt/IHC_cilia_RPParams.tc-1; a0=1; +b0=1+ a1; +% high pass (i.e. low pass reversed) +IHCciliaFilter_b=[a0 a1]; IHCciliaFilter_a=b0; +% figure(9), freqz(IHCciliaFilter_b, IHCciliaFilter_a) + +IHCciliaBndry=cell(nBFs,1); + +% IHC apical conductance (Boltzman function) +IHC_C= IHC_cilia_RPParams.C; +IHCu0= IHC_cilia_RPParams.u0; +IHCu1= IHC_cilia_RPParams.u1; +IHCs0= IHC_cilia_RPParams.s0; +IHCs1= IHC_cilia_RPParams.s1; +IHCGmax= IHC_cilia_RPParams.Gmax; +IHCGa= IHC_cilia_RPParams.Ga; % (leakage) + +IHCGu0 = IHCGa+IHCGmax./(1+exp(IHCu0/IHCs0).*(1+exp(IHCu1/IHCs1))); + +% Receptor potential +IHC_Cab= IHC_cilia_RPParams.Cab; +IHC_Gk= IHC_cilia_RPParams.Gk; +IHC_Et= IHC_cilia_RPParams.Et; +IHC_Ek= IHC_cilia_RPParams.Ek; +IHC_Ekp= IHC_Ek+IHC_Et*IHC_cilia_RPParams.Rpc; + +IHCrestingV= (IHC_Gk*IHC_Ekp+IHCGu0*IHC_Et)/(IHCGu0+IHC_Gk); + +IHC_Vnow= IHCrestingV*ones(nBFs,1); % initial voltage +IHC_RP= zeros(nBFs,segmentLength); + +% complete record of IHC receptor potential (V) +IHCciliaDisplacement= zeros(nBFs,segmentLength); +IHCoutput= zeros(nBFs,signalLength); +IHC_cilia_output= zeros(nBFs,signalLength); + + +%% pre-synapse --- +% Each BF is replicated using a different fiber type to make a 'channel' +% The number of channels is nBFs x nANfiberTypes +% Fiber types are specified in terms of tauCa +nANfiberTypes= length(IHCpreSynapseParams.tauCa); +tauCas= IHCpreSynapseParams.tauCa; +nChannels= nANfiberTypes*nBFs; +synapticCa= zeros(nChannels,segmentLength); + +% Calcium control (more calcium, greater release rate) +ECa=IHCpreSynapseParams.ECa; +gamma=IHCpreSynapseParams.gamma; +beta=IHCpreSynapseParams.beta; +tauM=IHCpreSynapseParams.tauM; +mICa=zeros(nChannels,segmentLength); +GmaxCa=IHCpreSynapseParams.GmaxCa; +synapse_z= IHCpreSynapseParams.z; +synapse_power=IHCpreSynapseParams.power; + +% tauCa vector is established across channels to allow vectorization +% (one tauCa per channel). Do not confuse with tauCas (one pre fiber type) +tauCa=repmat(tauCas, nBFs,1); +tauCa=reshape(tauCa, nChannels, 1); + +% presynapse startup values (vectors, length:nChannels) +% proportion (0 - 1) of Ca channels open at IHCrestingV +mICaCurrent=((1+beta^-1 * exp(-gamma*IHCrestingV))^-1)... + *ones(nBFs*nANfiberTypes,1); +% corresponding startup currents +ICaCurrent= (GmaxCa*mICaCurrent.^3) * (IHCrestingV-ECa); +CaCurrent= ICaCurrent.*tauCa; + +% vesicle release rate at startup (one per channel) +% kt0 is used only at initialisation +kt0= -synapse_z * CaCurrent.^synapse_power; + + +%% AN --- +% each row of the AN matrices represents one AN fiber +% The results computed either for probabiities *or* for spikes (not both) +% Spikes are necessary if CN and IC are to be computed +nFibersPerChannel= AN_IHCsynapseParams.numFibers; +nANfibers= nChannels*nFibersPerChannel; + +y=AN_IHCsynapseParams.y; +l=AN_IHCsynapseParams.l; +x=AN_IHCsynapseParams.x; +r=AN_IHCsynapseParams.r; +M=round(AN_IHCsynapseParams.M); + +% probability (NB initial 'P' on everything) +PAN_ydt = repmat(AN_IHCsynapseParams.y*dt, nChannels,1); +PAN_ldt = repmat(AN_IHCsynapseParams.l*dt, nChannels,1); +PAN_xdt = repmat(AN_IHCsynapseParams.x*dt, nChannels,1); +PAN_rdt = repmat(AN_IHCsynapseParams.r*dt, nChannels,1); +PAN_rdt_plus_ldt = PAN_rdt + PAN_ldt; +PAN_M=round(AN_IHCsynapseParams.M); + +% compute starting values +Pcleft = kt0* y* M ./ (y*(l+r)+ kt0* l); +Pavailable = Pcleft*(l+r)./kt0; +Preprocess = Pcleft*r/x; % canbe fractional + +ANprobability=zeros(nChannels,segmentLength); +ANprobRateOutput=zeros(nChannels,signalLength); +% special variables for monitoring synaptic cleft (specialists only) +savePavailableSeg=zeros(nChannels,segmentLength); +savePavailable=zeros(nChannels,signalLength); + +% spikes % ! ! ! ! ! ! ! ! +AN_refractory_period= AN_IHCsynapseParams.refractory_period; +lengthAbsRefractory= round(AN_refractory_period/ANdt); + +AN_ydt= repmat(AN_IHCsynapseParams.y*ANdt, nANfibers,1); +AN_ldt= repmat(AN_IHCsynapseParams.l*ANdt, nANfibers,1); +AN_xdt= repmat(AN_IHCsynapseParams.x*ANdt, nANfibers,1); +AN_rdt= repmat(AN_IHCsynapseParams.r*ANdt, nANfibers,1); +AN_rdt_plus_ldt= AN_rdt + AN_ldt; +AN_M= round(AN_IHCsynapseParams.M); + +% kt0 is initial release rate +% Establish as a vector (length=channel x number of fibers) +kt0= repmat(kt0', nFibersPerChannel, 1); +kt0=reshape(kt0, nANfibers,1); + +% starting values for reservoirs +AN_cleft = kt0* y* M ./ (y*(l+r)+ kt0* l); +AN_available = round(AN_cleft*(l+r)./kt0); %must be integer +AN_reprocess = AN_cleft*r/x; + +% output is in a logical array spikes = 1/ 0. +ANspikes= false(nANfibers,reducedSegmentLength); +ANoutput= false(nANfibers,reducedSignalLength); + + +%% CN (first brain stem nucleus - could be any subdivision of CN) +% Input to a CN neuorn is a random selection of AN fibers within a channel +% The number of AN fibers used is ANfibersFanInToCN +ANfibersFanInToCN=MacGregorMultiParams.fibersPerNeuron; +nCNneuronsPerChannel=MacGregorMultiParams.nNeuronsPerBF; +% CNtauGk (Potassium time constant) determines the rate of firing of +% the unit when driven hard by a DC input (not normally >350 sp/s) +CNtauGk=MacGregorMultiParams.tauGk; +ANavailableFibersPerChan=AN_IHCsynapseParams.numFibers; +nCNneurons=nCNneuronsPerChannel*nChannels; +% nCNneuronsPerFiberType= nCNneurons/nANfiberTypes; + +CNmembranePotential=zeros(nCNneurons,reducedSegmentLength); + +% establish which ANfibers (by name) feed into which CN nuerons +CNinputfiberLists=zeros(nChannels*nCNneuronsPerChannel, ANfibersFanInToCN); +unitNo=1; +for ch=1:nChannels + % Each channel contains a number of units =length(listOfFanInValues) + for idx=1:nCNneuronsPerChannel + fibersUsed=(ch-1)*ANavailableFibersPerChan + ... + ceil(rand(1,ANfibersFanInToCN)* ANavailableFibersPerChan); + CNinputfiberLists(unitNo,:)=fibersUsed; + unitNo=unitNo+1; + end +end + +% input to CN units +AN_PSTH=zeros(nCNneurons,reducedSegmentLength); + +% Generate CNalphaFunction function +% by which spikes are converted to post-synaptic currents +CNdendriteLPfreq= MacGregorMultiParams.dendriteLPfreq; +CNcurrentPerSpike=MacGregorMultiParams.currentPerSpike; +CNspikeToCurrentTau=1/(2*pi*CNdendriteLPfreq); +t=ANdt:ANdt:5*CNspikeToCurrentTau; +CNalphaFunction=... + (CNcurrentPerSpike/CNspikeToCurrentTau)*t.*exp(-t/CNspikeToCurrentTau); +% figure(98), plot(t,CNalphaFunction) +% working memory for implementing convolution +CNcurrentTemp=... + zeros(nCNneurons,reducedSegmentLength+length(CNalphaFunction)-1); +% trailing alphas are parts of humps carried forward to the next segment +CNtrailingAlphas=zeros(nCNneurons,length(CNalphaFunction)); + +CN_tauM=MacGregorMultiParams.tauM; +CN_tauTh=MacGregorMultiParams.tauTh; +CN_cap=MacGregorMultiParams.Cap; +CN_c=MacGregorMultiParams.c; +CN_b=MacGregorMultiParams.dGkSpike; +CN_Ek=MacGregorMultiParams.Ek; +CN_Eb= MacGregorMultiParams.Eb; +CN_Er=MacGregorMultiParams.Er; +CN_Th0= MacGregorMultiParams.Th0; +CN_E= zeros(nCNneurons,1); +CN_Gk= zeros(nCNneurons,1); +CN_Th= MacGregorMultiParams.Th0*ones(nCNneurons,1); +CN_Eb=CN_Eb.*ones(nCNneurons,1); +CN_Er=CN_Er.*ones(nCNneurons,1); +CNtimeSinceLastSpike=zeros(nCNneurons,1); +% tauGk is the main distinction between neurons +% in fact they are all the same in the standard model +tauGk=repmat(CNtauGk,nChannels*nCNneuronsPerChannel,1); + +CN_PSTH=zeros(nChannels,reducedSegmentLength); +CNoutput=false(nCNneurons,reducedSignalLength); + + +%% MacGregor (IC - second nucleus) -------- +nICcells=nChannels; % one cell per channel + +ICspikeWidth=0.00015; % this may need revisiting +epochsPerSpike=round(ICspikeWidth/ANdt); +if epochsPerSpike<1 + error(['MacGregorMulti: sample rate too low to support ' ... + num2str(ICspikeWidth*1e6) ' microsec spikes']); +end + +% short names +IC_tauM=MacGregorParams.tauM; +IC_tauGk=MacGregorParams.tauGk; +IC_tauTh=MacGregorParams.tauTh; +IC_cap=MacGregorParams.Cap; +IC_c=MacGregorParams.c; +IC_b=MacGregorParams.dGkSpike; +IC_Th0=MacGregorParams.Th0; +IC_Ek=MacGregorParams.Ek; +IC_Eb= MacGregorParams.Eb; +IC_Er=MacGregorParams.Er; + +IC_E=zeros(nICcells,1); +IC_Gk=zeros(nICcells,1); +IC_Th=IC_Th0*ones(nICcells,1); + +% Dendritic filtering, all spikes are replaced by CNalphaFunction functions +ICdendriteLPfreq= MacGregorParams.dendriteLPfreq; +ICcurrentPerSpike=MacGregorParams.currentPerSpike; +ICspikeToCurrentTau=1/(2*pi*ICdendriteLPfreq); +t=ANdt:ANdt:3*ICspikeToCurrentTau; +IC_CNalphaFunction= (ICcurrentPerSpike / ... + ICspikeToCurrentTau)*t.*exp(-t / ICspikeToCurrentTau); +% figure(98), plot(t,IC_CNalphaFunction) + +% working space for implementing alpha function +ICcurrentTemp=... + zeros(nICcells,reducedSegmentLength+length(IC_CNalphaFunction)-1); +ICtrailingAlphas=zeros(nICcells, length(IC_CNalphaFunction)); + +ICfiberTypeRates=zeros(nANfiberTypes,reducedSignalLength); +ICoutput=false(nChannels,reducedSignalLength); + +ICmembranePotential=zeros(nICcells,reducedSegmentLength); +ICmembraneOutput=zeros(nICcells,signalLength); + + +%% Main program %% %% %% %% %% %% %% %% %% %% %% %% %% %% + +% Compute the entire model for each segment +segmentStartPTR=1; +reducedSegmentPTR=1; % when sampling rate is reduced +while segmentStartPTR<signalLength + segmentEndPTR=segmentStartPTR+segmentLength-1; + % shorter segments after speed up. + shorterSegmentEndPTR=reducedSegmentPTR+reducedSegmentLength-1; + + iputPressureSegment=inputSignal... + (:,segmentStartPTR:segmentStartPTR+segmentLength-1); + + % segment debugging plots + % figure(98) + % plot(segmentTime,iputPressureSegment), title('signalSegment') + + + % OME ---------------------- + + % OME Stage 1: external resonances. Add to inputSignal pressure wave + y=iputPressureSegment; + for n=1:nOMEExtFilters + % any number of resonances can be used + [x OMEExtFilterBndry{n}] = ... + filter(ExtFilter_b{n},ExtFilter_a{n},... + iputPressureSegment, OMEExtFilterBndry{n}); + x= x* OMEgainScalars(n); + % This is a parallel resonance so add it + y=y+x; + end + iputPressureSegment=y; + OMEextEarPressure(segmentStartPTR:segmentEndPTR)= iputPressureSegment; + + % OME stage 2: convert input pressure (velocity) to + % tympanic membrane(TM) displacement using low pass filter + [TMdisplacementSegment OME_TMdisplacementBndry] = ... + filter(TMdisp_b,TMdisp_a,iputPressureSegment, ... + OME_TMdisplacementBndry); + % and save it + TMoutput(segmentStartPTR:segmentEndPTR)= TMdisplacementSegment; + + % OME stage 3: middle ear high pass effect to simulate stapes inertia + [stapesDisplacement OMEhighPassBndry] = ... + filter(stapesDisp_b,stapesDisp_a,TMdisplacementSegment, ... + OMEhighPassBndry); + + % OME stage 4: apply stapes scalar + stapesDisplacement=stapesDisplacement*stapesScalar; + + % OME stage 5: acoustic reflex stapes attenuation + % Attenuate the TM response using feedback from LSR fiber activity + if segmentStartPTR>efferentDelayPts + stapesDisplacement= stapesDisplacement.*... + ARattenuation(segmentStartPTR-efferentDelayPts:... + segmentEndPTR-efferentDelayPts); + end + + % segment debugging plots + % figure(98) + % plot(segmentTime, stapesDisplacement), title ('stapesDisplacement') + + % and save + OMEoutput(segmentStartPTR:segmentEndPTR)= stapesDisplacement; + + % needed for parallel processing + if segmentStartPTR>efferentDelayPts + MOCatt=MOCattenuation(:, segmentStartPTR-efferentDelayPts:... + segmentEndPTR-efferentDelayPts); + else % no MOC available yet + MOCatt=ones(nBFs, segmentLength); + end + %% BM ------------------------------ + % Each BM location is computed separately + parfor BFno=1:nBFs + + % *linear* path + % repeats used to avoid parallel processin problems + linOutput = stapesDisplacement * linGAIN; % linear gain + [linOutput GTlinBdry1out{BFno}] = ... + filter(GTlin_b(BFno,:), GTlin_a(BFno,:), linOutput, GTlinBdry1{BFno}); + [linOutput GTlinBdry2out{BFno}] = ... + filter(GTlin_b(BFno,:), GTlin_a(BFno,:), linOutput, GTlinBdry2{BFno}); + [linOutput GTlinBdry3out{BFno}] = ... + filter(GTlin_b(BFno,:), GTlin_a(BFno,:), linOutput, GTlinBdry3{BFno}); + + % *nonLinear* path + % efferent attenuation (0 <> 1) + MOC=MOCatt(BFno); + + + % first gammatone filter + [nonlinOutput firstGTnonlinBdry1out{BFno}] = ... + filter(GTnonlin_b(BFno,:), GTnonlin_a(BFno,:), ... + stapesDisplacement, firstGTnonlinBdry1{BFno}); + + [nonlinOutput firstGTnonlinBdry2out{BFno}] = ... + filter(GTnonlin_b(BFno,:), GTnonlin_a(BFno,:), ... + stapesDisplacement, firstGTnonlinBdry2{BFno}); + + [nonlinOutput firstGTnonlinBdry3out{BFno}] = ... + filter(GTnonlin_b(BFno,:), GTnonlin_a(BFno,:), ... + stapesDisplacement, firstGTnonlinBdry3{BFno}); + + % broken stick instantaneous compression + % nonlinear gain is weakend by MOC before applied to BM response + y= nonlinOutput.*(MOC* DRNLa); % linear section. + % compress those parts of the signal above the compression + % threshold + abs_x = abs(nonlinOutput); + idx=find(abs_x>DRNLcompressionThreshold); + if ~isempty(idx)>0 + y(idx)=sign(nonlinOutput(idx)).*... + (DRNLb*abs_x(idx).^DRNLc); + end + nonlinOutput=y; + + % second filter removes distortion products + [nonlinOutput secondGTnonlinBdry1out{BFno}] = ... + filter(GTnonlin_b(BFno,:), GTnonlin_a(BFno,:), ... + nonlinOutput, secondGTnonlinBdry1{BFno}); + + [nonlinOutput secondGTnonlinBdry2out{BFno}] = ... + filter(GTnonlin_b(BFno,:), GTnonlin_a(BFno,:), ... + nonlinOutput, secondGTnonlinBdry2{BFno}); + + [nonlinOutput secondGTnonlinBdry3out{BFno}] = ... + filter(GTnonlin_b(BFno,:), GTnonlin_a(BFno,:), ... + nonlinOutput, secondGTnonlinBdry3{BFno}); + + + % combine the two paths to give the DRNL displacement + DRNLresponse(BFno,:)=linOutput+nonlinOutput; + end % BF + GTlinBdry1=GTlinBdry1out; + GTlinBdry2=GTlinBdry2out; + GTlinBdry3=GTlinBdry3out; + firstGTnonlinBdry1=firstGTnonlinBdry1out; + firstGTnonlinBdry2=firstGTnonlinBdry2out; + firstGTnonlinBdry3=firstGTnonlinBdry3out; + secondGTnonlinBdry1=secondGTnonlinBdry1out; + secondGTnonlinBdry2=secondGTnonlinBdry2out; + secondGTnonlinBdry3=secondGTnonlinBdry3out; + % segment debugging plots + % figure(98) + % if size(DRNLresponse,1)>3 + % imagesc(DRNLresponse) % matrix display + % title('DRNLresponse'); % single or double channel response + % else + % plot(segmentTime, DRNLresponse) + % end + + % and save it + DRNLoutput(:, segmentStartPTR:segmentEndPTR)= DRNLresponse; + + + %% IHC ------------------------------------ + % BM displacement to IHCciliaDisplacement is a high-pass filter + % because of viscous coupling + for idx=1:nBFs + [IHCciliaDisplacement(idx,:) IHCciliaBndry{idx}] = ... + filter(IHCciliaFilter_b,IHCciliaFilter_a, ... + DRNLresponse(idx,:), IHCciliaBndry{idx}); + end + + % apply scalar + IHCciliaDisplacement=IHCciliaDisplacement* IHC_C; + + % and save it + IHC_cilia_output(:,segmentStartPTR:segmentStartPTR+segmentLength-1)=... + IHCciliaDisplacement; + + % compute apical conductance + G=IHCGmax./(1+exp(-(IHCciliaDisplacement-IHCu0)/IHCs0).*... + (1+exp(-(IHCciliaDisplacement-IHCu1)/IHCs1))); + Gu=G + IHCGa; + + % Compute receptor potential + for idx=1:segmentLength + IHC_Vnow=IHC_Vnow+ (-Gu(:, idx).*(IHC_Vnow-IHC_Et)-... + IHC_Gk*(IHC_Vnow-IHC_Ekp))* dt/IHC_Cab; + IHC_RP(:,idx)=IHC_Vnow; + end + + % segment debugging plots + % if size(IHC_RP,1)>3 + % surf(IHC_RP), shading interp, title('IHC_RP') + % else + % plot(segmentTime, IHC_RP) + % end + + % and save it + IHCoutput(:, segmentStartPTR:segmentStartPTR+segmentLength-1)=IHC_RP; + + + %% synapse ----------------------------- + % Compute the vesicle release rate for each fiber type at each BF + % replicate IHC_RP for each fiber type + Vsynapse=repmat(IHC_RP, nANfiberTypes,1); + + % look-up table of target fraction channels open for a given IHC_RP + mICaINF= 1./( 1 + exp(-gamma * Vsynapse) /beta); + % fraction of channel open - apply time constant + for idx=1:segmentLength + % mICaINF is the current 'target' value of mICa + mICaCurrent=mICaCurrent+(mICaINF(:,idx)-mICaCurrent)*dt./tauM; + mICa(:,idx)=mICaCurrent; + end + + ICa= (GmaxCa* mICa.^3) .* (Vsynapse- ECa); + + for idx=1:segmentLength + CaCurrent=CaCurrent + ICa(:,idx)*dt - CaCurrent*dt./tauCa; + synapticCa(:,idx)=CaCurrent; + end + synapticCa=-synapticCa; % treat IHCpreSynapseParams as positive substance + + % NB vesicleReleaseRate is /s and is independent of dt + vesicleReleaseRate = synapse_z * synapticCa.^synapse_power; % rate + + % segment debugging plots + % if size(vesicleReleaseRate,1)>3 + % surf(vesicleReleaseRate), shading interp, title('vesicleReleaseRate') + % else + % plot(segmentTime, vesicleReleaseRate) + % end + + + %% AN + switch AN_spikesOrProbability + case 'probability' + % No refractory effect is applied + for t = 1:segmentLength; + M_Pq=PAN_M-Pavailable; + M_Pq(M_Pq<0)=0; + Preplenish = M_Pq .* PAN_ydt; + Pejected = Pavailable.* vesicleReleaseRate(:,t)*dt; + Preprocessed = M_Pq.*Preprocess.* PAN_xdt; + + ANprobability(:,t)= min(Pejected,1); + reuptakeandlost= PAN_rdt_plus_ldt .* Pcleft; + reuptake= PAN_rdt.* Pcleft; + + Pavailable= Pavailable+ Preplenish- Pejected+ Preprocessed; + Pcleft= Pcleft + Pejected - reuptakeandlost; + Preprocess= Preprocess + reuptake - Preprocessed; + Pavailable(Pavailable<0)=0; + savePavailableSeg(:,t)=Pavailable; % synapse tracking + end + % and save it as *rate* + ANrate=ANprobability/dt; + ANprobRateOutput(:, segmentStartPTR:... + segmentStartPTR+segmentLength-1)= ANrate; + % monitor synapse contents (only sometimes used) + savePavailable(:, segmentStartPTR:segmentStartPTR+segmentLength-1)=... + savePavailableSeg; + + % Estimate efferent effects. ARattenuation (0 <> 1) + % acoustic reflex + ARAttSeg=mean(ANrate(1:nBFs,:),1); %LSR channels are 1:nBF + % smooth + [ARAttSeg, ARboundaryProb] = ... + filter(ARfilt_b, ARfilt_a, ARAttSeg, ARboundaryProb); + ARAttSeg=ARAttSeg-ARrateThreshold; + ARAttSeg(ARAttSeg<0)=0; % prevent negative strengths + ARattenuation(segmentStartPTR:segmentEndPTR)=... + (1-ARrateToAttenuationFactorProb.* ARAttSeg); + + % MOC attenuation + % within-channel HSR response only + HSRbegins=nBFs*(nANfiberTypes-1)+1; + rates=ANrate(HSRbegins:end,:); + for idx=1:nBFs + [smoothedRates, MOCprobBoundary{idx}] = ... + filter(MOCfilt_b, MOCfilt_a, rates(idx,:), ... + MOCprobBoundary{idx}); + smoothedRates=smoothedRates-MOCrateThreshold; + smoothedRates(smoothedRates<0)=0; + MOCattenuation(idx,segmentStartPTR:segmentEndPTR)= ... + (1- smoothedRates* rateToAttenuationFactorProb); + end + MOCattenuation(MOCattenuation<0)=0.001; + + + case 'spikes' + ANtimeCount=0; + % implement speed upt + for t = ANspeedUpFactor:ANspeedUpFactor:segmentLength; + ANtimeCount=ANtimeCount+1; + % convert release rate to probabilities + releaseProb=vesicleReleaseRate(:,t)*ANdt; + % releaseProb is the release probability per channel + % but each channel has many synapses + releaseProb=repmat(releaseProb',nFibersPerChannel,1); + releaseProb=reshape(releaseProb, nFibersPerChannel*nChannels,1); + + % AN_available=round(AN_available); % vesicles must be integer, (?needed) + M_q=AN_M- AN_available; % number of missing vesicles + M_q(M_q<0)= 0; % cannot be less than 0 + + % AN_N1 converts probability to discrete events + % it considers each event that might occur + % (how many vesicles might be released) + % and returns a count of how many were released + + % slow line +% probabilities= 1-(1-releaseProb).^AN_available; + probabilities= 1-intpow((1-releaseProb), AN_available); + ejected= probabilities> rand(length(AN_available),1); + + reuptakeandlost = AN_rdt_plus_ldt .* AN_cleft; + reuptake = AN_rdt.* AN_cleft; + + % slow line +% probabilities= 1-(1-AN_reprocess.*AN_xdt).^M_q; + probabilities= 1-intpow((1-AN_reprocess.*AN_xdt), M_q); + reprocessed= probabilities>rand(length(M_q),1); + + % slow line +% probabilities= 1-(1-AN_ydt).^M_q; + probabilities= 1-intpow((1-AN_ydt), M_q); + + replenish= probabilities>rand(length(M_q),1); + + AN_available = AN_available + replenish - ejected ... + + reprocessed; + AN_cleft = AN_cleft + ejected - reuptakeandlost; + AN_reprocess = AN_reprocess + reuptake - reprocessed; + + % ANspikes is logical record of vesicle release events>0 + ANspikes(:, ANtimeCount)= ejected; + end % t + + % zero any events that are preceded by release events ... + % within the refractory period + % The refractory period consist of two periods + % 1) the absolute period where no spikes occur + % 2) a relative period where a spike may occur. This relative + % period is realised as a variable length interval + % where the length is chosen at random + % (esentially a linear ramp up) + + % Andreas has a fix for this + for t = 1:ANtimeCount-2*lengthAbsRefractory; + % identify all spikes across fiber array at time (t) + % idx is a list of channels where spikes occurred + % ?? try sparse matrices? + idx=find(ANspikes(:,t)); + for j=idx % consider each spike + % specify variable refractory period + % between abs and 2*abs refractory period + nPointsRefractory=lengthAbsRefractory+... + round(rand*lengthAbsRefractory); + % disable spike potential for refractory period + % set all values in this range to 0 + ANspikes(j,t+1:t+nPointsRefractory)=0; + end + end %t + + % segment debugging + % plotInstructions.figureNo=98; + % plotInstructions.displaydt=ANdt; + % plotInstructions.numPlots=1; + % plotInstructions.subPlotNo=1; + % UTIL_plotMatrix(ANspikes, plotInstructions); + + % and save it. NB, AN is now on 'speedUp' time + ANoutput(:, reducedSegmentPTR: shorterSegmentEndPTR)=ANspikes; + + + %% CN Macgregor first neucleus ------------------------------- + % input is from AN so ANdt is used throughout + % Each CNneuron has a unique set of input fibers selected + % at random from the available AN fibers (CNinputfiberLists) + + % Create the dendritic current for that neuron + % First get input spikes to this neuron + synapseNo=1; + for ch=1:nChannels + for idx=1:nCNneuronsPerChannel + % determine candidate fibers for this unit + fibersUsed=CNinputfiberLists(synapseNo,:); + % ANpsth has a bin width of dt + % (just a simple sum across fibers) + AN_PSTH(synapseNo,:) = ... + sum(ANspikes(fibersUsed,:), 1); + synapseNo=synapseNo+1; + end + end + + % One alpha function per spike + [alphaRows alphaCols]=size(CNtrailingAlphas); + + for unitNo=1:nCNneurons + CNcurrentTemp(unitNo,:)= ... + conv(AN_PSTH(unitNo,:),CNalphaFunction); + end + % add post-synaptic current left over from previous segment + CNcurrentTemp(:,1:alphaCols)=... + CNcurrentTemp(:,1:alphaCols)+ CNtrailingAlphas; + + % take post-synaptic current for this segment + CNcurrentInput= CNcurrentTemp(:, 1:reducedSegmentLength); + + % trailingalphas are the ends of the alpha functions that + % spill over into the next segment + CNtrailingAlphas= ... + CNcurrentTemp(:, reducedSegmentLength+1:end); + + if CN_c>0 + % variable threshold condition (slow) + for t=1:reducedSegmentLength + CNtimeSinceLastSpike=CNtimeSinceLastSpike-dts; + s=CN_E>CN_Th & CNtimeSinceLastSpike<0 ; + CNtimeSinceLastSpike(s)=0.0005; % 0.5 ms for sodium spike + dE =(-CN_E/CN_tauM + ... + CNcurrentInput(:,t)/CN_cap+(CN_Gk/CN_cap).*(CN_Ek-CN_E))*dt; + dGk=-CN_Gk*dt./tauGk + CN_b*s; + dTh=-(CN_Th-CN_Th0)*dt/CN_tauTh + CN_c*s; + CN_E=CN_E+dE; + CN_Gk=CN_Gk+dGk; + CN_Th=CN_Th+dTh; + CNmembranePotential(:,t)=CN_E+s.*(CN_Eb-CN_E)+CN_Er; + end + else + % static threshold (faster) + for t=1:reducedSegmentLength + CNtimeSinceLastSpike=CNtimeSinceLastSpike-dt; + s=CN_E>CN_Th0 & CNtimeSinceLastSpike<0 ; % =1 if both conditions met + CNtimeSinceLastSpike(s)=0.0005; % 0.5 ms for sodium spike + dE = (-CN_E/CN_tauM + ... + CNcurrentInput(:,t)/CN_cap+(CN_Gk/CN_cap).*(CN_Ek-CN_E))*dt; + dGk=-CN_Gk*dt./tauGk +CN_b*s; + CN_E=CN_E+dE; + CN_Gk=CN_Gk+dGk; + % add spike to CN_E and add resting potential (-60 mV) + CNmembranePotential(:,t)=CN_E+s.*(CN_Eb-CN_E)+CN_Er; + end + end + + % extract spikes. A spike is a substantial upswing in voltage + CN_spikes=CNmembranePotential> -0.01; + + % now remove any spike that is immediately followed by a spike + % NB 'find' works on columns (whence the transposing) + CN_spikes=CN_spikes'; + idx=find(CN_spikes); + idx=idx(1:end-1); + CN_spikes(idx+1)=0; + CN_spikes=CN_spikes'; + + % segment debugging + % plotInstructions.figureNo=98; + % plotInstructions.displaydt=ANdt; + % plotInstructions.numPlots=1; + % plotInstructions.subPlotNo=1; + % UTIL_plotMatrix(CN_spikes, plotInstructions); + + % and save it + CNoutput(:, reducedSegmentPTR:shorterSegmentEndPTR)=... + CN_spikes; + + + %% IC ---------------------------------------------- + % MacGregor or some other second order neurons + + % combine CN neurons in same channel, i.e. same BF & same tauCa + % to generate inputs to single IC unit + channelNo=0; + for idx=1:nCNneuronsPerChannel:nCNneurons-nCNneuronsPerChannel+1; + channelNo=channelNo+1; + CN_PSTH(channelNo,:)=... + sum(CN_spikes(idx:idx+nCNneuronsPerChannel-1,:)); + end + + [alphaRows alphaCols]=size(ICtrailingAlphas); + for ICneuronNo=1:nICcells + ICcurrentTemp(ICneuronNo,:)= ... + conv(CN_PSTH(ICneuronNo,:), IC_CNalphaFunction); + end + + % add the unused current from the previous convolution + ICcurrentTemp(:,1:alphaCols)=ICcurrentTemp(:,1:alphaCols)... + + ICtrailingAlphas; + % take what is required and keep the trailing part for next time + inputCurrent=ICcurrentTemp(:, 1:reducedSegmentLength); + ICtrailingAlphas=ICcurrentTemp(:, reducedSegmentLength+1:end); + + if IC_c==0 + % faster computation when threshold is stable (C==0) + for t=1:reducedSegmentLength + s=IC_E>IC_Th0; + dE = (-IC_E/IC_tauM + inputCurrent(:,t)/IC_cap +... + (IC_Gk/IC_cap).*(IC_Ek-IC_E))*dt; + dGk=-IC_Gk*dt/IC_tauGk +IC_b*s; + IC_E=IC_E+dE; + IC_Gk=IC_Gk+dGk; + ICmembranePotential(:,t)=IC_E+s.*(IC_Eb-IC_E)+IC_Er; + end + else + % threshold is changing (IC_c>0; e.g. bushy cell) + for t=1:reducedSegmentLength + dE = (-IC_E/IC_tauM + ... + inputCurrent(:,t)/IC_cap + (IC_Gk/IC_cap)... + .*(IC_Ek-IC_E))*dt; + IC_E=IC_E+dE; + s=IC_E>IC_Th; + ICmembranePotential(:,t)=IC_E+s.*(IC_Eb-IC_E)+IC_Er; + dGk=-IC_Gk*dt/IC_tauGk +IC_b*s; + IC_Gk=IC_Gk+dGk; + + % After a spike, the threshold is raised + % otherwise it settles to its baseline + dTh=-(IC_Th-Th0)*dt/IC_tauTh +s*IC_c; + IC_Th=IC_Th+dTh; + end + end + + ICspikes=ICmembranePotential> -0.01; + % now remove any spike that is immediately followed by a spike + % NB 'find' works on columns (whence the transposing) + ICspikes=ICspikes'; + idx=find(ICspikes); + idx=idx(1:end-1); + ICspikes(idx+1)=0; + ICspikes=ICspikes'; + + nCellsPerTau= nICcells/nANfiberTypes; + firstCell=1; + lastCell=nCellsPerTau; + for tauCount=1:nANfiberTypes + % separate rates according to fiber types + ICfiberTypeRates(tauCount, ... + reducedSegmentPTR:shorterSegmentEndPTR)=... + sum(ICspikes(firstCell:lastCell, :))... + /(nCellsPerTau*ANdt); + firstCell=firstCell+nCellsPerTau; + lastCell=lastCell+nCellsPerTau; + end + ICoutput(:, reducedSegmentPTR:shorterSegmentEndPTR)=ICspikes; + + if nBFs==1 % single channel + x= repmat(ICmembranePotential(1,:), ANspeedUpFactor,1); + x= reshape(x,1,segmentLength); + if nANfiberTypes>1 % save HSR and LSR + y= repmat(ICmembranePotential(end,:), ANspeedUpFactor,1); + y= reshape(y,1,segmentLength); + x=[x; y]; + end + ICmembraneOutput(:, segmentStartPTR:segmentEndPTR)= x; + end + + % estimate efferent effects. + % ARis based on LSR units. LSR channels are 1:nBF + if nANfiberTypes>1 % AR is multi-channel only + ARAttSeg=sum(ICspikes(1:nBFs,:),1)/ANdt; + [ARAttSeg, ARboundary] = ... + filter(ARfilt_b, ARfilt_a, ARAttSeg, ARboundary); + ARAttSeg=ARAttSeg-ARrateThreshold; + ARAttSeg(ARAttSeg<0)=0; % prevent negative strengths + % scale up to dt from ANdt + x= repmat(ARAttSeg, ANspeedUpFactor,1); + x=reshape(x,1,segmentLength); + ARattenuation(segmentStartPTR:segmentEndPTR)=... + (1-ARrateToAttenuationFactor* x); + ARattenuation(ARattenuation<0)=0.001; + else + % single channel model; disable AR + ARattenuation(segmentStartPTR:segmentEndPTR)=... + ones(1,segmentLength); + end + + % MOC attenuation using HSR response only + % Separate MOC effect for each BF + HSRbegins=nBFs*(nANfiberTypes-1)+1; + rates=ICspikes(HSRbegins:end,:)/ANdt; + for idx=1:nBFs + [smoothedRates, MOCboundary{idx}] = ... + filter(MOCfilt_b, MOCfilt_a, rates(idx,:), ... + MOCboundary{idx}); + MOCattSegment(idx,:)=smoothedRates; + % expand timescale back to model dt from ANdt + x= repmat(MOCattSegment(idx,:), ANspeedUpFactor,1); + x= reshape(x,1,segmentLength); + MOCattenuation(idx,segmentStartPTR:segmentEndPTR)= ... + (1- MOCrateToAttenuationFactor* x); + end + MOCattenuation(MOCattenuation<0)=0.04; + % segment debugging + % plotInstructions.figureNo=98; + % plotInstructions.displaydt=ANdt; + % plotInstructions.numPlots=1; + % plotInstructions.subPlotNo=1; + % UTIL_plotMatrix(ICspikes, plotInstructions); + + end % AN_spikesOrProbability + segmentStartPTR=segmentStartPTR+segmentLength; + reducedSegmentPTR=reducedSegmentPTR+reducedSegmentLength; + + +end % segment + +path(restorePath)
--- a/multithreshold 1.46/expGUI_MT.m Fri May 27 17:50:40 2011 +0100 +++ b/multithreshold 1.46/expGUI_MT.m Tue May 31 09:13:07 2011 +0100 @@ -46,7 +46,7 @@ % Edit the above text to modify the response to help expGUI_MT -% Last Modified by GUIDE v2.5 07-Apr-2011 07:20:00 +% Last Modified by GUIDE v2.5 29-May-2011 16:02:02 % Begin initialization code - DO NOT EDIT gui_Singleton = 1; @@ -314,7 +314,7 @@ paradigm=paradigmNames{chosenOption}; experiment.paradigm=paradigm; -%Paradigm: read in all relevant parameters + %Paradigm: read in all relevant parameters % a file must exist with this name 'paradigm_<paradigm>' % 'handles' are only occasionally used addpath ('paradigms') @@ -327,6 +327,16 @@ end rmpath ('paradigms') +switch experiment.paradigm + % identify masker free paradigms + case {'training', 'discomfort','absThreshold', 'absThreshold_8',... + 'absThreshold_16','TENtest', 'threshold_duration','SRT'... + 'paradigm_thr_IFMC'} + experiment.maskerInUse=0; + otherwise + experiment.maskerInUse=1; +end + % if a variable is subject to change, specify list of values here % eg. stimulusParameters.targetFrequency=betweenRuns.variableList1; cmd=['stimulusParameters.' ... @@ -423,7 +433,7 @@ aShowRelevantObjects(handles) -% % backgroundType +% % backgroundType popup idx= find(strcmp(stimulusParameters.backgroundType, backgroundTypes)); set(handles.popupmenuBackgroundType,'value', idx) set(handles.editBackgroundLevel,'string', ... @@ -432,6 +442,7 @@ % ---------------------------------------------- aShowRelevantObjects function aShowRelevantObjects(handles) global experiment stimulusParameters +% called from aShowRelevantObjects % always on set(handles.edittargetLevel, 'visible', 'on') set(handles.edittargetDuration, 'visible', 'on') @@ -442,15 +453,7 @@ set(handles.editStatsModel, 'visible', 'on') set(handles.textStatsModel, 'visible', 'on') set(handles.pushbuttonStop, 'visible', 'on') - set(handles.pushbuttonOME, 'visible', 'off') - set(handles.pushbuttonBM, 'visible', 'off') - set(handles.pushbuttonRP, 'visible', 'off') - set(handles.pushbuttonAN, 'visible', 'off') - set(handles.pushbuttonRF, 'visible', 'off') - set(handles.pushbuttonSYN, 'visible', 'off') - set(handles.pushbuttonFM, 'visible', 'off') - set(handles.pushbuttonParams, 'visible', 'off') - set(handles.pushbuttonSingleShot, 'visible', 'off') +showModelPushButtons(handles, 0) set(handles.editCatchTrialRate, 'visible', 'off') set(handles.textCatchTrials, 'visible', 'off') set(handles.editcalibrationdB, 'visible', 'off') @@ -459,44 +462,25 @@ set(handles.textCue, 'visible', 'off') set(handles.editMusicLevel,'visible', 'off') set(handles.textMusicLevel,'visible', 'off') - set(handles.editMAPplot,'visible', 'off') - set(handles.textMAPplot,'visible', 'off') case {'MAPmodel', 'MAPmodelListen', 'MAPmodelMultiCh', 'MAPmodelSingleCh'} set(handles.popupmenuCueNoCue, 'visible', 'off') set(handles.editStatsModel, 'visible', 'off') set(handles.textStatsModel, 'visible', 'off') set(handles.pushbuttonStop, 'visible', 'on') - set(handles.pushbuttonOME, 'visible', 'on') - set(handles.pushbuttonBM, 'visible', 'on') - set(handles.pushbuttonRP, 'visible', 'on') - set(handles.pushbuttonAN, 'visible', 'on') - set(handles.pushbuttonRF, 'visible', 'on') - set(handles.pushbuttonSYN, 'visible', 'on') - set(handles.pushbuttonFM, 'visible', 'on') - set(handles.pushbuttonParams, 'visible', 'on') - set(handles.pushbuttonSingleShot, 'visible', 'on') +showModelPushButtons(handles, 1) set(handles.editcalibrationdB, 'visible', 'off') set(handles.textcalibration, 'visible', 'off') set(handles.textCue, 'visible', 'off') set(handles.editMusicLevel,'visible', 'off') set(handles.textMusicLevel,'visible', 'off') - set(handles.editMAPplot,'visible', 'on') - set(handles.textMAPplot,'visible', 'on') otherwise + % i.e. using real subjects (left, right, diotic, dichotic) set(handles.editStatsModel, 'visible', 'off') set(handles.textStatsModel, 'visible', 'off') set(handles.pushbuttonStop, 'visible', 'off') - set(handles.pushbuttonOME, 'visible', 'off') - set(handles.pushbuttonBM, 'visible', 'off') - set(handles.pushbuttonRP, 'visible', 'off') - set(handles.pushbuttonAN, 'visible', 'off') - set(handles.pushbuttonRF, 'visible', 'off') - set(handles.pushbuttonSYN, 'visible', 'off') - set(handles.pushbuttonFM, 'visible', 'off') - set(handles.pushbuttonParams, 'visible', 'off') - set(handles.pushbuttonSingleShot, 'visible', 'off') +showModelPushButtons(handles, 0) set(handles.editCatchTrialRate, 'visible', 'on') set(handles.textCatchTrials, 'visible', 'on') set(handles.editcalibrationdB, 'visible', 'on') @@ -505,8 +489,6 @@ set(handles.textCue, 'visible', 'on') set(handles.editMusicLevel,'visible', 'on') set(handles.textMusicLevel,'visible', 'on') - set(handles.editMAPplot,'visible', 'off') - set(handles.textMAPplot,'visible', 'off') end switch experiment.threshEstMethod @@ -533,48 +515,42 @@ set(handles.textcueTestDifference,'visible', 'off') end +% show/ remove masker related boxes +if ~experiment.maskerInUse + set(handles.editmaskerDuration,'visible', 'off') + set(handles.editmaskerLevel,'visible', 'off') + set(handles.editmaskerRelativeFrequency,'visible', 'off') + set(handles.editgapDuration,'visible', 'off') + set(handles.textmaskerDuration,'visible', 'off') + set(handles.textmaskerLevel,'visible', 'off') + set(handles.textmaskerRelativeFrequency,'visible', 'off') + set(handles.textgapDuration,'visible', 'off') + set(handles.popupmenuMaskerType,'visible', 'off') + set(handles.textMaskerType,'visible', 'off') -% show/ remove masker related boxes -switch experiment.paradigm - % masker free paradigms - case {'training', 'discomfort','absThreshold', 'absThreshold_8',... - 'TENtest', 'threshold_duration','SRT'} - set(handles.editmaskerDuration,'visible', 'off') - set(handles.editmaskerLevel,'visible', 'off') - set(handles.editmaskerRelativeFrequency,'visible', 'off') - set(handles.editgapDuration,'visible', 'off') - set(handles.textmaskerDuration,'visible', 'off') - set(handles.textmaskerLevel,'visible', 'off') - set(handles.textmaskerRelativeFrequency,'visible', 'off') - set(handles.textgapDuration,'visible', 'off') - set(handles.popupmenuMaskerType,'visible', 'off') - set(handles.textMaskerType,'visible', 'off') - - % paradigms with maskers - case {'forwardMasking','forwardMaskingD','trainingIFMC', 'TMC','TMC_16ms', ... - 'TMCmodel','IFMC','IFMC_8ms', 'IFMC_16ms','GOM','overShoot','overShootB',... - 'gapDetection', 'psychometric', 'FMreProbe'} - set(handles.editmaskerDuration,'visible', 'on') - set(handles.editmaskerLevel,'visible', 'on') - set(handles.editmaskerRelativeFrequency,'visible', 'on') - set(handles.editgapDuration,'visible', 'on') - set(handles.popupmenuMaskerType,'visible', 'on') - set(handles.textmaskerDuration,'visible', 'on') - set(handles.textmaskerLevel,'visible', 'on') - set(handles.textmaskerRelativeFrequency,'visible', 'on') - set(handles.textgapDuration,'visible', 'on') - set(handles.popupmenuMaskerType,'visible', 'on') - set(handles.textMaskerType,'visible', 'on') - % masker relative frequency /type - chosenOption=get(handles.popupmenuMaskerType,'value'); - maskerTypes=get(handles.popupmenuMaskerType,'string'); - maskerType=maskerTypes{chosenOption}; - switch maskerType - case 'tone' - set(handles.editmaskerRelativeFrequency,'visible', 'on') - otherwise - set(handles.editmaskerRelativeFrequency,'visible', 'off') - end + % paradigms with maskers +else + set(handles.editmaskerDuration,'visible', 'on') + set(handles.editmaskerLevel,'visible', 'on') + set(handles.editmaskerRelativeFrequency,'visible', 'on') + set(handles.editgapDuration,'visible', 'on') + set(handles.popupmenuMaskerType,'visible', 'on') + set(handles.textmaskerDuration,'visible', 'on') + set(handles.textmaskerLevel,'visible', 'on') + set(handles.textmaskerRelativeFrequency,'visible', 'on') + set(handles.textgapDuration,'visible', 'on') + set(handles.popupmenuMaskerType,'visible', 'on') + set(handles.textMaskerType,'visible', 'on') + % masker relative frequency /type + chosenOption=get(handles.popupmenuMaskerType,'value'); + maskerTypes=get(handles.popupmenuMaskerType,'string'); + maskerType=maskerTypes{chosenOption}; + switch maskerType + case 'tone' + set(handles.editmaskerRelativeFrequency,'visible', 'on') + otherwise + set(handles.editmaskerRelativeFrequency,'visible', 'off') + end end eval(['set(handles.edit' stimulusParameters.WRVname ... @@ -601,16 +577,94 @@ set(handles.textBGlevel,'visible','on') end +% ------------------------------------------------ showModelPushButtons +function showModelPushButtons(handles,on) +if on + set(handles.pushbuttonOME, 'visible', 'on') + set(handles.pushbuttonBM, 'visible', 'on') + set(handles.pushbuttonRP, 'visible', 'on') + set(handles.pushbuttonAN, 'visible', 'on') + set(handles.pushbuttonPhLk, 'visible', 'on') + set(handles.pushbuttonSYN, 'visible', 'on') + set(handles.pushbuttonFM, 'visible', 'on') + set(handles.pushbuttonParams, 'visible', 'on') + set(handles.pushbuttonSingleShot, 'visible', 'on') + set(handles.editMAPplot,'visible', 'on') + set(handles.textMAPplot,'visible', 'on') +else + set(handles.pushbuttonOME, 'visible', 'off') + set(handles.pushbuttonBM, 'visible', 'off') + set(handles.pushbuttonRP, 'visible', 'off') + set(handles.pushbuttonAN, 'visible', 'off') + set(handles.pushbuttonPhLk, 'visible', 'off') + set(handles.pushbuttonSYN, 'visible', 'off') + set(handles.pushbuttonFM, 'visible', 'off') + set(handles.pushbuttonParams, 'visible', 'off') + set(handles.pushbuttonSingleShot, 'visible', 'off') + set(handles.editMAPplot,'visible', 'off') + set(handles.textMAPplot,'visible', 'off') + +end % ------------------------------------------------ pushbuttonRun_Callback function pushbuttonRun_Callback(hObject, eventdata, handles) global checkForPreviousGUI % holds screen positioning across repeated calls -global experiment +global experiment betweenRuns paradigmNames checkForPreviousGUI.GUIposition=get(handles.figure1,'position'); experiment.singleShot=0; -run (handles) -experiment.stop=0; +switch experiment.paradigm + case 'thr_IFMC' + %% special option for two successive and linked measurements + experiment.paradigm='threshold_16ms'; + set(handles.editstopCriteriaBox,'string','30') % nTrials + run (handles) + + % use these threshold for IFMC + thresholds16ms=betweenRuns.thresholds; + optionNo=strmatch('IFMC_16ms',paradigmNames); + set(handles.popupmenuParadigm,'value',optionNo); + aParadigmSelection(handles) + set(handles.edittargetLevel,'string', thresholds16ms+10); + set(handles.editstopCriteriaBox,'string','30') % nTrials + pause(.1) + run (handles) + + % reset original paradigm + optionNo=strmatch('thr_IFMC',paradigmNames); + set(handles.popupmenuParadigm,'value',optionNo); + aParadigmSelection(handles) + + case 'thr_TMC' + %% special option for two successive and linked measurements + experiment.paradigm='threshold_16ms'; + set(handles.edittargetDuration,'string', num2str(0.25)) + set(handles.editstopCriteriaBox,'string','30') % nTrials + run (handles) + + set(handles.edittargetDuration,'string', num2str(0.016)) + set(handles.editstopCriteriaBox,'string','30') % nTrials + run (handles) + + % use these threshold for TMC + thresholds16ms=betweenRuns.thresholds; + optionNo=strmatch('TMC_16ms',paradigmNames); + set(handles.popupmenuParadigm,'value',optionNo); + aParadigmSelection(handles) + set(handles.edittargetLevel,'string', thresholds16ms+10); + set(handles.editstopCriteriaBox,'string','30') % nTrials + pause(.1) + run (handles) + + % reset original paradigm + optionNo=strmatch('thr_TMC',paradigmNames); + set(handles.popupmenuParadigm,'value',optionNo); + aParadigmSelection(handles) + + otherwise + run (handles) + experiment.stop=0; +end function run (handles) global experiment expGUIhandles stimulusParameters @@ -926,10 +980,10 @@ return end +% set the target level in advance for every wnticipated run switch experiment.paradigm - % case { 'TMC', 'TMCmodel','IFMC'} - case {'trainingIFMC', 'TMC','TMC_16ms', 'TMC - ELP', 'IFMC', 'IFMC_16ms'} - + case {'trainingIFMC', 'TMC','TMC_16ms', 'TMC - ELP', ... + 'IFMC', 'IFMC_16ms'} % For TMC and IFMC multiple target levels can be set if length(stimulusParameters.targetLevel)==1 betweenRuns.targetLevels=... @@ -938,6 +992,7 @@ *length(betweenRuns.variableList2)); elseif length(stimulusParameters.targetLevel)==... length(betweenRuns.variableList2) + % only one level specified, so use this throughout x=stimulusParameters.targetLevel; x=repmat(x,length(betweenRuns.variableList1),1); x=reshape(x,1,length(betweenRuns.variableList1)*length(betweenRuns.variableList2)); @@ -1309,15 +1364,18 @@ testBM(stimulusParameters.targetFrequency, experiment.name); function pushbuttonAN_Callback(hObject, eventdata, handles) +global stimulusParameters aReadAndCheckParameterBoxes(handles); % now carry out tests showPSTHs=0; -testAN +targetFrequency=stimulusParameters.targetFrequency(1); +BFlist=targetFrequency; -function pushbuttonRF_Callback(hObject, eventdata, handles) +testAN(targetFrequency,BFlist); + +function pushbuttonPhLk_Callback(hObject, eventdata, handles) aReadAndCheckParameterBoxes(handles); -showPSTHs=1; -testRF +testPhaseLocking function pushbuttonSYN_Callback(hObject, eventdata, handles) aReadAndCheckParameterBoxes(handles);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/multithreshold 1.46/hs_err_pid6852.log Tue May 31 09:13:07 2011 +0100 @@ -0,0 +1,288 @@ +# +# An unexpected error has been detected by Java Runtime Environment: +# +# EXCEPTION_ACCESS_VIOLATION (0xc0000005) at pc=0x0000000011cc5110, pid=6852, tid=5540 +# +# Java VM: Java HotSpot(TM) 64-Bit Server VM (1.6.0-b105 mixed mode) +# Problematic frame: +# C [awt.dll+0x185110] +# +# If you would like to submit a bug report, please visit: +# http://java.sun.com/webapps/bugreport/crash.jsp +# + +--------------- T H R E A D --------------- + +Current thread (0x0000000011841800): JavaThread "AWT-EventQueue-0" [_thread_in_native, id=5540] + +siginfo: ExceptionCode=0xc0000005, reading address 0xffffffffffffffff + +Registers: +EAX=0x800000b26e4dd767, EBX=0x0000000000000001, ECX=0x0000000030776480, EDX=0x0000000011a012d0 +ESP=0x000000003126e600, EBP=0x000000001183a360, ESI=0x0000000011841990, EDI=0x0000000000000000 +EIP=0x0000000011cc5110, EFLAGS=0x0000000000010206 + +Top of Stack: (sp=0x000000003126e600) +0x000000003126e600: 0000000011841800 0000000000000000 +0x000000003126e610: 000000003126e750 0000000015e00818 +0x000000003126e620: 0000000000000001 0000000011cc6dd3 +0x000000003126e630: 0000000000000001 000000001183a360 +0x000000003126e640: 0000000011841990 0000000000000000 +0x000000003126e650: 0000000011a012d0 00000000169f4028 +0x000000003126e660: 0000000000000004 00000000152cd810 +0x000000003126e670: 0000000000000000 0000000000000000 +0x000000003126e680: 0000000000000001 0000000000000102 +0x000000003126e690: 00000000122a308e 0000000000000001 +0x000000003126e6a0: 000000003126e728 00000000122afc24 +0x000000003126e6b0: 0000000000000000 0000000026630e70 +0x000000003126e6c0: 0000000000000000 0000000015e00818 +0x000000003126e6d0: 0000000000000001 0000000000000000 +0x000000003126e6e0: 000000003126e6d0 000000003126e6e8 +0x000000003126e6f0: 0000000000000000 000000003126e750 + +Instructions: (pc=0x0000000011cc5110) +0x0000000011cc5100: e8 5b 72 f9 ff 48 8b 0d 04 9a 08 00 48 8b 04 d8 +0x0000000011cc5110: 48 8b 58 30 e8 67 72 f9 ff 48 8b c3 48 83 c4 20 + + +Stack: [0x00000000311f0000,0x0000000031270000), sp=0x000000003126e600, free space=505k +Native frames: (J=compiled Java code, j=interpreted, Vv=VM code, C=native code) +C [awt.dll+0x185110] + +Java frames: (J=compiled Java code, j=interpreted, Vv=VM code) +j sun.awt.Win32GraphicsConfig.getBounds(I)Ljava/awt/Rectangle;+0 +j sun.awt.Win32GraphicsConfig.getBounds()Ljava/awt/Rectangle;+8 +j com.mathworks.mwswing.WindowUtils.getVirtualScreenBounds()Ljava/awt/Rectangle;+98 +j com.mathworks.mwswing.MJUtilities.getVirtualScreenBounds()Ljava/awt/Rectangle;+0 +j com.mathworks.mwswing.desk.DTSingleClientFrame.refineLocation(Lcom/mathworks/mwswing/desk/DTLocation;)Lcom/mathworks/mwswing/desk/DTFloatingLocation;+495 +j com.mathworks.mwswing.desk.DTSingleClientFrame.addClient(Lcom/mathworks/mwswing/desk/DTClient;Lcom/mathworks/mwswing/desk/DTLocation;)V+297 +j com.mathworks.mwswing.desk.Desktop.setClientShowing(Lcom/mathworks/mwswing/desk/DTClient;ZLcom/mathworks/mwswing/desk/DTLocation;ZZ)V+956 +j com.mathworks.mwswing.desk.Desktop.setClientShowing(Lcom/mathworks/mwswing/desk/DTClient;ZLcom/mathworks/mwswing/desk/DTLocation;Z)V+7 +j com.mathworks.mwswing.desk.Desktop.setClientDocked(Lcom/mathworks/mwswing/desk/DTClient;Z)V+296 +j com.mathworks.mwswing.desk.DTClient$ClientUndockAction.actionPerformed(Ljava/awt/event/ActionEvent;)V+12 +j javax.swing.AbstractButton.fireActionPerformed(Ljava/awt/event/ActionEvent;)V+84 +j javax.swing.AbstractButton$Handler.actionPerformed(Ljava/awt/event/ActionEvent;)V+5 +j javax.swing.DefaultButtonModel.fireActionPerformed(Ljava/awt/event/ActionEvent;)V+35 +j javax.swing.DefaultButtonModel.setPressed(Z)V+117 +j javax.swing.plaf.basic.BasicButtonListener.mouseReleased(Ljava/awt/event/MouseEvent;)V+35 +j java.awt.AWTEventMulticaster.mouseReleased(Ljava/awt/event/MouseEvent;)V+8 +j java.awt.AWTEventMulticaster.mouseReleased(Ljava/awt/event/MouseEvent;)V+8 +j java.awt.Component.processMouseEvent(Ljava/awt/event/MouseEvent;)V+64 +j javax.swing.JComponent.processMouseEvent(Ljava/awt/event/MouseEvent;)V+23 +j com.mathworks.mwswing.MJButton.processMouseEvent(Ljava/awt/event/MouseEvent;)V+37 +j java.awt.Component.processEvent(Ljava/awt/AWTEvent;)V+81 +j java.awt.Container.processEvent(Ljava/awt/AWTEvent;)V+18 +J java.awt.Component.dispatchEventImpl(Ljava/awt/AWTEvent;)V +J java.awt.Container.dispatchEventImpl(Ljava/awt/AWTEvent;)V +J java.awt.LightweightDispatcher.retargetMouseEvent(Ljava/awt/Component;ILjava/awt/event/MouseEvent;)V +J java.awt.Container.dispatchEventImpl(Ljava/awt/AWTEvent;)V +j java.awt.Window.dispatchEventImpl(Ljava/awt/AWTEvent;)V+19 +J java.awt.EventDispatchThread.pumpOneEventForFilters(I)Z +J java.awt.EventDispatchThread.pumpEventsForFilter(ILjava/awt/Conditional;Ljava/awt/EventFilter;)V +j java.awt.EventDispatchThread.pumpEventsForHierarchy(ILjava/awt/Conditional;Ljava/awt/Component;)V+11 +j java.awt.EventDispatchThread.pumpEvents(ILjava/awt/Conditional;)V+4 +j java.awt.EventDispatchThread.pumpEvents(Ljava/awt/Conditional;)V+3 +j java.awt.EventDispatchThread.run()V+9 +v ~StubRoutines::call_stub + +--------------- P R O C E S S --------------- + +Java Threads: ( => current thread ) + 0x0000000031567800 JavaThread "Timer-7" [_thread_blocked, id=4716] + 0x0000000031568000 JavaThread "Image Fetcher 0" daemon [_thread_blocked, id=6432] + 0x0000000031566c00 JavaThread "Inactive RequestProcessor thread [Was:TimedSoftReference/org.openide.util.TimedSoftReference]" daemon [_thread_blocked, id=7912] + 0x0000000031565000 JavaThread "Thread-237" [_thread_blocked, id=5244] + 0x0000000031568800 JavaThread "Thread-230" [_thread_blocked, id=5112] + 0x0000000031567400 JavaThread "Prefs Updater" [_thread_blocked, id=6764] + 0x0000000031566400 JavaThread "Thread-10" [_thread_blocked, id=7204] + 0x0000000031565c00 JavaThread "Thread-9" [_thread_blocked, id=6880] + 0x0000000031565800 JavaThread "Thread-8" [_thread_blocked, id=6964] + 0x0000000031564800 JavaThread "Active Reference Queue Daemon" daemon [_thread_blocked, id=7260] + 0x0000000031564000 JavaThread "Timer-3" daemon [_thread_blocked, id=4892] + 0x000000003155d000 JavaThread "Timer-2" daemon [_thread_blocked, id=6948] + 0x0000000030c13000 JavaThread "TimerQueue" daemon [_thread_blocked, id=5596] +=>0x0000000011841800 JavaThread "AWT-EventQueue-0" [_thread_in_native, id=5540] + 0x00000000307b0800 JavaThread "AWT-Windows" daemon [_thread_in_native, id=5832] + 0x00000000307afc00 JavaThread "AWT-Shutdown" [_thread_blocked, id=4548] + 0x00000000307af400 JavaThread "Java2D Disposer" daemon [_thread_blocked, id=5040] + 0x0000000011779000 JavaThread "Timer-0" [_thread_blocked, id=6884] + 0x00000000111f1c00 JavaThread "Low Memory Detector" daemon [_thread_blocked, id=4740] + 0x00000000111ef800 JavaThread "CompilerThread1" daemon [_thread_blocked, id=7472] + 0x00000000111e7c00 JavaThread "CompilerThread0" daemon [_thread_blocked, id=6800] + 0x00000000111e6800 JavaThread "Attach Listener" daemon [_thread_blocked, id=7688] + 0x00000000111c7400 JavaThread "Finalizer" daemon [_thread_blocked, id=7524] + 0x0000000003def400 JavaThread "Reference Handler" daemon [_thread_blocked, id=2228] + 0x0000000003d2fc00 JavaThread "main" [_thread_in_native, id=5124] + +Other Threads: + 0x0000000003dee400 VMThread [id=7268] + 0x00000000111f7400 WatcherThread [id=7276] + +VM state:not at safepoint (normal execution) + +VM Mutex/Monitor currently owned by a thread: None + +Heap + PSYoungGen total 11392K, used 691K [0x00000000265a0000, 0x00000000273a0000, 0x00000000296a0000) + eden space 10624K, 6% used [0x00000000265a0000,0x000000002664cc98,0x0000000027000000) + from space 768K, 0% used [0x0000000027040000,0x0000000027040000,0x0000000027100000) + to space 1856K, 0% used [0x00000000271d0000,0x00000000271d0000,0x00000000273a0000) + PSOldGen total 50944K, used 35225K [0x000000001d2a0000, 0x0000000020460000, 0x00000000265a0000) + object space 50944K, 69% used [0x000000001d2a0000,0x000000001f5066a8,0x0000000020460000) + PSPermGen total 58624K, used 44603K [0x00000000152a0000, 0x0000000018be0000, 0x000000001d2a0000) + object space 58624K, 76% used [0x00000000152a0000,0x0000000017e2ef18,0x0000000018be0000) + +Dynamic libraries: +0x0000000140000000 - 0x0000000140138000 C:\Program Files\MATLAB\R2008a\bin\win64\MATLAB.exe +0x0000000077a60000 - 0x0000000077c09000 C:\Windows\SYSTEM32\ntdll.dll +0x0000000077940000 - 0x0000000077a5f000 C:\Windows\system32\kernel32.dll +0x000007fefdbd0000 - 0x000007fefdc3b000 C:\Windows\system32\KERNELBASE.dll +0x0000000180000000 - 0x0000000180401000 C:\Program Files\MATLAB\R2008a\bin\win64\libut.dll +0x000007fefe1b0000 - 0x000007fefe1c7000 C:\Windows\system32\imagehlp.dll +0x000007fefdf80000 - 0x000007fefe01f000 C:\Windows\system32\msvcrt.dll +0x0000000077c30000 - 0x0000000077c37000 C:\Windows\system32\PSAPI.DLL +0x00000000010d0000 - 0x00000000010f8000 C:\Program Files\MATLAB\R2008a\bin\win64\LIBEXPAT.dll +0x0000000073fa0000 - 0x0000000074069000 C:\Windows\WinSxS\amd64_microsoft.vc80.crt_1fc8b3b9a1e18e3b_8.0.50727.5592_none_88e45feb2faab9ce\MSVCR80.dll +0x000000004a800000 - 0x000000004a917000 C:\Program Files\MATLAB\R2008a\bin\win64\icuuc36.dll +0x000007feffc30000 - 0x000007feffd0b000 C:\Windows\system32\ADVAPI32.dll +0x000007feff4e0000 - 0x000007feff4ff000 C:\Windows\SYSTEM32\sechost.dll +0x000007feff140000 - 0x000007feff26d000 C:\Windows\system32\RPCRT4.dll +0x0000000001120000 - 0x0000000001123000 C:\Program Files\MATLAB\R2008a\bin\win64\icudt36.dll +0x000000004ab00000 - 0x000000004ab0f000 C:\Program Files\MATLAB\R2008a\bin\win64\icuio36.dll +0x00000000012b0000 - 0x00000000013b6000 C:\Program Files\MATLAB\R2008a\bin\win64\icuin36.dll +0x0000000072920000 - 0x0000000072a29000 C:\Windows\WinSxS\amd64_microsoft.vc80.crt_1fc8b3b9a1e18e3b_8.0.50727.5592_none_88e45feb2faab9ce\MSVCP80.dll +0x0000000077840000 - 0x000000007793a000 C:\Windows\system32\USER32.dll +0x000007feff5a0000 - 0x000007feff607000 C:\Windows\system32\GDI32.dll +0x000007fefe1a0000 - 0x000007fefe1ae000 C:\Windows\system32\LPK.dll +0x000007feffb60000 - 0x000007feffc29000 C:\Windows\system32\USP10.dll +0x00000000013c0000 - 0x0000000001534000 C:\Program Files\MATLAB\R2008a\bin\win64\libmwservices.dll +0x0000000001540000 - 0x00000000015b0000 C:\Program Files\MATLAB\R2008a\bin\win64\libmx.dll +0x0000000001180000 - 0x0000000001197000 C:\Program Files\MATLAB\R2008a\bin\win64\zlib1.dll +0x00000000015b0000 - 0x0000000001658000 C:\Program Files\MATLAB\R2008a\bin\win64\libmwmathutil.dll +0x0000000001660000 - 0x00000000016b5000 C:\Program Files\MATLAB\R2008a\bin\win64\mpath.dll +0x00000000016d0000 - 0x00000000016f1000 C:\Program Files\MATLAB\R2008a\bin\win64\mlutil.dll +0x000007fef9630000 - 0x000007fef96d0000 C:\Windows\WinSxS\amd64_microsoft.windows.common-controls_6595b64144ccf1df_5.82.7601.17514_none_a4d6a923711520a9\COMCTL32.dll +0x000007fefdd80000 - 0x000007fefde17000 C:\Windows\system32\comdlg32.dll +0x000007feff660000 - 0x000007feff6d1000 C:\Windows\system32\SHLWAPI.dll +0x000007fefe1d0000 - 0x000007fefef58000 C:\Windows\system32\SHELL32.dll +0x000007fefce80000 - 0x000007fefce96000 C:\Windows\system32\NETAPI32.dll +0x000007fefce70000 - 0x000007fefce7c000 C:\Windows\system32\netutils.dll +0x000007fefd500000 - 0x000007fefd523000 C:\Windows\system32\srvcli.dll +0x000007fefce50000 - 0x000007fefce65000 C:\Windows\system32\wkscli.dll +0x000007feff610000 - 0x000007feff65d000 C:\Windows\system32\WS2_32.dll +0x000007feff6e0000 - 0x000007feff6e8000 C:\Windows\system32\NSI.dll +0x0000000001710000 - 0x0000000001765000 C:\Program Files\MATLAB\R2008a\bin\win64\mcr.dll +0x0000000001780000 - 0x00000000017a5000 C:\Program Files\MATLAB\R2008a\bin\win64\iqm.dll +0x00000000017c0000 - 0x00000000017e1000 C:\Program Files\MATLAB\R2008a\bin\win64\bridge.dll +0x0000000001800000 - 0x0000000001811000 C:\Program Files\MATLAB\R2008a\bin\win64\libmex.dll +0x0000000001830000 - 0x00000000018bc000 C:\Program Files\MATLAB\R2008a\bin\win64\m_dispatcher.dll +0x00000000018d0000 - 0x00000000018f5000 C:\Program Files\MATLAB\R2008a\bin\win64\datasvcs.dll +0x0000000012000000 - 0x0000000012295000 C:\Program Files\MATLAB\R2008a\bin\win64\xerces-c_2_7.dll +0x0000000001920000 - 0x00000000021b1000 C:\Program Files\MATLAB\R2008a\bin\win64\m_interpreter.dll +0x00000000021d0000 - 0x0000000002201000 C:\Program Files\MATLAB\R2008a\bin\win64\libmat.dll +0x0000000002220000 - 0x0000000002325000 C:\Program Files\MATLAB\R2008a\bin\win64\libhdf5.dll +0x0000000002330000 - 0x000000000239f000 C:\Program Files\MATLAB\R2008a\bin\win64\profiler.dll +0x00000000023b0000 - 0x00000000023ba000 C:\Program Files\MATLAB\R2008a\bin\win64\libmwmathrng.dll +0x00000000023d0000 - 0x00000000023ea000 C:\Program Files\MATLAB\R2008a\bin\win64\m_pcodeio.dll +0x0000000002400000 - 0x000000000244a000 C:\Program Files\MATLAB\R2008a\bin\win64\m_ir.dll +0x0000000002460000 - 0x0000000002a1b000 C:\Program Files\MATLAB\R2008a\bin\win64\m_parser.dll +0x0000000002a30000 - 0x0000000002a42000 C:\Program Files\MATLAB\R2008a\bin\win64\ir_xfmr.dll +0x0000000002a60000 - 0x0000000002c7b000 C:\Program Files\MATLAB\R2008a\bin\win64\mcos.dll +0x0000000002c90000 - 0x0000000002c9c000 C:\Program Files\MATLAB\R2008a\bin\win64\mtok.dll +0x0000000002cb0000 - 0x0000000002cd0000 C:\Program Files\MATLAB\R2008a\bin\win64\m_pcodegen.dll +0x000007fefad70000 - 0x000007fefae95000 C:\Windows\system32\dbghelp.dll +0x0000000002ce0000 - 0x0000000002cf0000 C:\Program Files\MATLAB\R2008a\bin\win64\boost_thread-vc80-mt-1_34_1.dll +0x0000000002d00000 - 0x0000000002dc0000 C:\Program Files\MATLAB\R2008a\bin\win64\udd.dll +0x0000000002dd0000 - 0x0000000002f12000 C:\Program Files\MATLAB\R2008a\bin\win64\libmwgui.dll +0x0000000002f30000 - 0x000000000316a000 C:\Program Files\MATLAB\R2008a\bin\win64\hg.dll +0x0000000003180000 - 0x00000000031d6000 C:\Program Files\MATLAB\R2008a\bin\win64\jmi.dll +0x00000000031f0000 - 0x000000000322e000 C:\Program Files\MATLAB\R2008a\bin\win64\libmwhardcopy.dll +0x0000000003240000 - 0x000000000329c000 C:\Program Files\MATLAB\R2008a\bin\win64\libuij.dll +0x00000000032b0000 - 0x000000000353c000 C:\Program Files\MATLAB\R2008a\bin\win64\numerics.dll +0x0000000003550000 - 0x000000000355c000 C:\Program Files\MATLAB\R2008a\bin\win64\libmwblas.dll +0x0000000003570000 - 0x000000000357f000 C:\Program Files\MATLAB\R2008a\bin\win64\libmwbinder.dll +0x0000000003590000 - 0x00000000035b4000 C:\Program Files\MATLAB\R2008a\bin\win64\libmwlapack.dll +0x00000000035d0000 - 0x00000000035db000 C:\Program Files\MATLAB\R2008a\bin\win64\libmwfftw.dll +0x00000000035f0000 - 0x0000000003625000 C:\Program Files\MATLAB\R2008a\bin\win64\libmwrookfastbp.dll +0x0000000003640000 - 0x000000000366e000 C:\Program Files\MATLAB\R2008a\bin\win64\libmwma57.dll +0x0000000010000000 - 0x00000000100d3000 C:\Program Files\MATLAB\R2008a\bin\win64\libifcoremd.dll +0x0000000003680000 - 0x000000000389d000 C:\Program Files\MATLAB\R2008a\bin\win64\libmmd.dll +0x00000000038a0000 - 0x00000000038a9000 C:\Program Files\MATLAB\R2008a\bin\win64\libmwcsparse.dll +0x00000000038c0000 - 0x000000000398a000 C:\Program Files\MATLAB\R2008a\bin\win64\libmwumfpack.dll +0x00000000039a0000 - 0x00000000039ad000 C:\Program Files\MATLAB\R2008a\bin\win64\libmwamd.dll +0x00000000039c0000 - 0x0000000003a52000 C:\Program Files\MATLAB\R2008a\bin\win64\libmwcholmod.dll +0x0000000003a70000 - 0x0000000003a7c000 C:\Program Files\MATLAB\R2008a\bin\win64\libmwcolamd.dll +0x0000000003a90000 - 0x0000000003b49000 C:\Program Files\MATLAB\R2008a\bin\win64\uiw.dll +0x0000000003b60000 - 0x0000000003b6a000 C:\Program Files\MATLAB\R2008a\bin\win64\uinone.dll +0x0000000071490000 - 0x000000007162c000 C:\Windows\WinSxS\amd64_microsoft.vc80.mfc_1fc8b3b9a1e18e3b_8.0.50727.5592_none_8448f49f328da8c3\MFC80.DLL +0x000007fefbcd0000 - 0x000007fefbd41000 C:\Windows\system32\WINSPOOL.DRV +0x000007feff6f0000 - 0x000007feff8f3000 C:\Windows\system32\ole32.dll +0x000007feff400000 - 0x000007feff4d7000 C:\Windows\system32\OLEAUT32.dll +0x0000000003b80000 - 0x0000000003c10000 C:\Program Files\MATLAB\R2008a\bin\win64\udd_mi.dll +0x0000000003c20000 - 0x0000000003c38000 C:\Program Files\MATLAB\R2008a\bin\win64\mwoles05.DLL +0x0000000003c50000 - 0x0000000003cb9000 C:\Program Files\MATLAB\R2008a\bin\win64\comcli.dll +0x0000000072900000 - 0x0000000072920000 C:\Windows\WinSxS\amd64_microsoft.vc80.atl_1fc8b3b9a1e18e3b_8.0.50727.5592_none_8a1e1b372ed7b012\ATL80.DLL +0x0000000003cd0000 - 0x0000000003cde000 C:\Program Files\MATLAB\R2008a\bin\win64\mlautoregister.dll +0x000007fefde20000 - 0x000007fefde4e000 C:\Windows\system32\IMM32.DLL +0x000007feff2f0000 - 0x000007feff3f9000 C:\Windows\system32\MSCTF.dll +0x000000006fa00000 - 0x000000006fa3f000 C:\PROGRA~2\Sophos\SOPHOS~1\SOPHOS~2.DLL +0x0000000007240000 - 0x0000000007c0e000 C:\Program Files\MATLAB\R2008a\bin\win64\mkl.dll +0x0000000005810000 - 0x000000000585b000 C:\Program Files\MATLAB\R2008a\bin\win64\libguide40.dll +0x0000000005860000 - 0x0000000005868000 C:\Program Files\MATLAB\R2008a\bin\win64\mklcompat.dll +0x0000000007c10000 - 0x00000000081b7000 C:\Program Files\MATLAB\R2008a\bin\win64\mllapack.dll +0x0000000006c40000 - 0x0000000006d34000 C:\Program Files\MATLAB\R2008a\bin\win64\libfftw3i.dll +0x0000000006d40000 - 0x0000000006e2e000 C:\Program Files\MATLAB\R2008a\bin\win64\libfftw3f.dll +0x000007fefd9b0000 - 0x000007fefd9bf000 C:\Windows\system32\profapi.dll +0x000007fefd8e0000 - 0x000007fefd8ef000 C:\Windows\system32\CRYPTBASE.dll +0x000007fefc4d0000 - 0x000007fefc6c4000 C:\Windows\WinSxS\amd64_microsoft.windows.common-controls_6595b64144ccf1df_6.0.7601.17514_none_fa396087175ac9ac\comctl32.dll +0x000007fefef60000 - 0x000007feff137000 C:\Windows\system32\SETUPAPI.dll +0x000007fefdca0000 - 0x000007fefdcd6000 C:\Windows\system32\CFGMGR32.dll +0x000007fefdc80000 - 0x000007fefdc9a000 C:\Windows\system32\DEVOBJ.dll +0x000007feff500000 - 0x000007feff599000 C:\Windows\system32\CLBCatQ.DLL +0x000007fefbf30000 - 0x000007fefc05c000 C:\Windows\system32\propsys.dll +0x000007fefb4e0000 - 0x000007fefb50d000 C:\Windows\system32\ntmarta.dll +0x000007feffd10000 - 0x000007feffd62000 C:\Windows\system32\WLDAP32.dll +0x000007fefa3d0000 - 0x000007fefa3f7000 C:\Windows\system32\iphlpapi.dll +0x000007fefa3c0000 - 0x000007fefa3cb000 C:\Windows\system32\WINNSI.DLL +0x000007fefd0a0000 - 0x000007fefd0fb000 C:\Windows\system32\DNSAPI.dll +0x000007fef9e00000 - 0x000007fef9e11000 C:\Windows\system32\dhcpcsvc6.DLL +0x000007fef9db0000 - 0x000007fef9dc8000 C:\Windows\system32\dhcpcsvc.DLL +0x000007fefd850000 - 0x000007fefd875000 C:\Windows\system32\SspiCli.dll +0x0000000074fe0000 - 0x0000000074fe3000 C:\Windows\system32\icmp.Dll +0x000000000f820000 - 0x000000000fd71000 C:\Program Files\MATLAB\R2008a\sys\java\jre\win64\jre1.6.0\bin\server\jvm.dll +0x000007fefb5c0000 - 0x000007fefb5fb000 C:\Windows\system32\WINMM.dll +0x0000000006aa0000 - 0x0000000006aaa000 C:\Program Files\MATLAB\R2008a\sys\java\jre\win64\jre1.6.0\bin\hpi.dll +0x00000000082c0000 - 0x00000000082ce000 C:\Program Files\MATLAB\R2008a\sys\java\jre\win64\jre1.6.0\bin\verify.dll +0x00000000082d0000 - 0x00000000082f7000 C:\Program Files\MATLAB\R2008a\sys\java\jre\win64\jre1.6.0\bin\java.dll +0x0000000008300000 - 0x0000000008312000 C:\Program Files\MATLAB\R2008a\sys\java\jre\win64\jre1.6.0\bin\zip.dll +0x000000000ff40000 - 0x000000000ffa5000 C:\Program Files\WIDCOMM\Bluetooth Software\btmmhook.dll +0x000000000ffb0000 - 0x000000000ffc6000 C:\Program Files\MATLAB\R2008a\bin\win64\nativejava.dll +0x000000000ffe0000 - 0x000000000fff6000 C:\Program Files\MATLAB\R2008a\bin\win64\nativejmi.dll +0x00000000112f0000 - 0x00000000112f7000 C:\Program Files\MATLAB\R2008a\bin\win64\nativeservices.dll +0x0000000011b40000 - 0x0000000011d90000 C:\Program Files\MATLAB\R2008a\sys\java\jre\win64\jre1.6.0\bin\awt.dll +0x0000000011f10000 - 0x0000000011f79000 C:\Program Files\MATLAB\R2008a\sys\java\jre\win64\jre1.6.0\bin\fontmanager.dll + +VM Arguments: +jvm_args: -Xss512k -XX:PermSize=32M -Xms64m -XX:NewRatio=3 -XX:MaxPermSize=128M -Xmx196m -XX:MaxDirectMemorySize=2147400000 -Dsun.java2d.noddraw=true -Dsun.awt.nopixfmt=true -Xshare:off -Xrs -Djava.library.path=C:\Program Files\MATLAB\R2008a\bin\win64 vfprintf abort +java_command: <unknown> +Launcher Type: generic + +Environment Variables: +CLASSPATH=.;C:\Program Files (x86)\Java\jre6\lib\ext\QTJava.zip +PATH=C:\Program Files (x86)\Nokia\PC Connectivity Solution\;C:\Windows\system32;C:\Windows;C:\Windows\System32\Wbem;C:\Windows\System32\WindowsPowerShell\v1.0\;C:\Program Files\Intel\WiFi\bin\;C:\Program Files\Common Files\Intel\WirelessCommon\;C:\Program Files\Intel\DMIX;C:\Program Files (x86)\NTRU Cryptosystems\NTRU TCG Software Stack\bin\;C:\Program Files\NTRU Cryptosystems\NTRU TCG Software Stack\bin\;C:\Program Files\Wave Systems Corp\Gemalto\Access Client\v5\;c:\Program Files\WIDCOMM\Bluetooth Software\;c:\Program Files\WIDCOMM\Bluetooth Software\syswow64;C:\Program Files (x86)\Common Files\Roxio Shared\DLLShared\;C:\Program Files (x86)\Common Files\Roxio Shared\10.0\DLLShared\;C:\Program Files (x86)\Common Files\Adobe\AGL;C:\Program Files\MATLAB\R2010b\bin;C:\Program Files\MATLAB\R2010a\bin;C:\Program Files\MATLAB\R2008a\bin;C:\Program Files\MATLAB\R2008a\bin\win64;C:\Program Files (x86)\QuickTime\QTSystem\;C:\Program Files\Microsoft Windows Performance Toolkit\ +USERNAME=rmeddis +OS=Windows_NT +PROCESSOR_IDENTIFIER=Intel64 Family 6 Model 37 Stepping 2, GenuineIntel + + + +--------------- S Y S T E M --------------- + +OS: Windows NT 6.1 Build 7601 Service Pack 1 + +CPU:total 4 em64t ht + +Memory: 4k page, physical 8181592k(5194804k free), swap 16361336k(13377632k free) + +vm_info: Java HotSpot(TM) 64-Bit Server VM (1.6.0-b105) for windows-amd64, built on Nov 29 2006 00:38:01 by "java_re" with unknown MS VC++:1400 +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/multithreshold 1.46/old files/MAPmodel.m Tue May 31 09:13:07 2011 +0100 @@ -0,0 +1,136 @@ +function [modelResponse, MacGregorResponse]=MAPmodel( MAPplot, method) + +global experiment stimulusParameters audio withinRuns +global outerMiddleEarParams DRNLParams AN_IHCsynapseParams + +savePath=path; +addpath('..\MAP') +modelResponse=[]; +MacGregorResponse=[]; + +% mono only (column vector) +audio=audio(:,1)'; + +% if stop button pressed earlier +if experiment.stop, return, end + +% -------------------------------------------------------------- run Model +MAPparamsName=experiment.name; +showPlotsAndDetails=experiment.MAPplot; +AN_spikesOrProbability='spikes'; + +% [response, method]=MAPsequenceSeg(audio, method, 1:8); +global ICoutput ANdt + MAP1_14(audio, 1/method.dt, method.nonlinCF,... + MAPparamsName, AN_spikesOrProbability); + +if showPlotsAndDetails + options.showModelParameters=0; + options.showModelOutput=1; + options.printFiringRates=1; + options.showACF=0; + options.showEfferent=1; + showMAP(options) +end + +% No response, probably caused by hitting 'stop' button +if isempty(ICoutput), return, end + +% MacGregor response is the sum total of all final stage spiking +MacGregorResponse= sum(ICoutput,1); % use IC + +% ---------------------------------------------------------- end model run + +dt=ANdt; +time=dt:dt:dt*length(MacGregorResponse); + +% group delay on unit response +MacGonsetDelay= 0.004; +MacGoffsetDelay= 0.022; + +% now find the response of the MacGregor model during the target presentation + group delay +switch experiment.threshEstMethod + case {'2I2AFC++', '2I2AFC+++'} + idx= time>stimulusParameters.testTargetBegins+MacGonsetDelay ... + & time<stimulusParameters.testTargetEnds+MacGoffsetDelay; + nSpikesTrueWindow=sum(MacGregorResponse(:,idx)); + idx=find(time>stimulusParameters.testNonTargetBegins+MacGonsetDelay ... + & time<stimulusParameters.testNonTargetEnds+MacGoffsetDelay); + nSpikesFalseWindow=sum(MacGregorResponse(:,idx)); + % nSpikesDuringTarget is +ve when more spikes are found + % in the target window + difference= nSpikesTrueWindow-nSpikesFalseWindow; + + if difference>0 + % hit + nSpikesDuringTarget=experiment.MacGThreshold+1; + elseif difference<0 + % miss (wrong choice) + nSpikesDuringTarget=experiment.MacGThreshold-1; + else + if rand>0.5 + % hit (random choice) + nSpikesDuringTarget=experiment.MacGThreshold+1; + else + % miss (random choice) + nSpikesDuringTarget=experiment.MacGThreshold-1; + end + end + disp(['level target dummy decision: ' ... + num2str([withinRuns.variableValue nSpikesTrueWindow ... + nSpikesFalseWindow nSpikesDuringTarget], '%4.0f') ] ) + + otherwise + % idx=find(time>stimulusParameters.testTargetBegins+MacGonsetDelay ... + % & time<stimulusParameters.testTargetEnds+MacGoffsetDelay); + % no delay at onset + idx=find(time>stimulusParameters.testTargetBegins +MacGonsetDelay... + & time<stimulusParameters.testTargetEnds+MacGoffsetDelay); + nSpikesDuringTarget=sum(MacGregorResponse(:,idx)); + + % find(MacGregorResponse)*dt-stimulusParameters.stimulusDelay + timeX=time(idx); +end + +% now find the response of the MacGregor model at the end of the masker +idx2=find(time>stimulusParameters.testTargetBegins-0.02 ... + & time<stimulusParameters.testTargetBegins); +if ~isempty(idx2) + maskerRate=mean(mean(MacGregorResponse(idx2))); +else + %e.g. no masker + maskerRate=0; +end + +if experiment.MAPplot + % add vertical lines to indicate target region + figure(99), subplot(6,1,6) + hold on + yL=get(gca,'YLim'); + plot([stimulusParameters.testTargetBegins + MacGonsetDelay ... + stimulusParameters.testTargetBegins + MacGonsetDelay],yL,'r') + plot([stimulusParameters.testTargetEnds + MacGoffsetDelay ... + stimulusParameters.testTargetEnds + MacGoffsetDelay],yL,'r') +end + +% specify unambiguous response +switch experiment.paradigm + case 'gapDetection' + gapResponse=(maskerRate-nSpikesDuringTarget)/maskerRate; + if gapResponse>0.2 + modelResponse=2; % gap detected + else + modelResponse=1; % gap not detected + end + [nSpikesDuringTarget maskerRate gapResponse modelResponse] + figure(22), plot(timeX,earObject(idx)) + otherwise + if nSpikesDuringTarget>experiment.MacGThreshold + modelResponse=2; % stimulus detected + else + modelResponse=1; % nothing heard (default) + end +end + + +path(savePath) \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/multithreshold 1.46/paradigms/paradigm_thr_IFMC.m Tue May 31 09:13:07 2011 +0100 @@ -0,0 +1,36 @@ +function paradigm_thr_IFMC(handles) +global stimulusParameters experiment betweenRuns + +paradigm_training(handles) % default + +betweenRuns.variableName1='targetFrequency'; +betweenRuns.variableList1=1000; +betweenRuns.variableList1=str2num(get(handles.edittargetFrequency,'string')); +betweenRuns.variableName2='targetDuration'; +betweenRuns.variableList2=0.016; +betweenRuns.randomizeSequence=1; % 'random sequence' + +% delay > masker > gap > target + + +stimulusParameters.targetType='tone'; +stimulusParameters.targetPhase='sin'; +stimulusParameters.targetFrequency=1000; +stimulusParameters.targetDuration=0.016; +stimulusParameters.targetLevel=stimulusParameters.WRVstartValues(1); + +stimulusParameters.rampDuration=0.004; + +experiment.stopCriteria2IFC=[75 3 5]; +experiment.singleIntervalMaxTrials=[20]; + + +% forced choice window interval +stimulusParameters.AFCsilenceDuration=0.5; + +% instructions to user +% single interval up/down no cue +stimulusParameters.instructions{1}= [{'YES if you hear the tone clearly'}, { }, { 'NO if not (or you are uncertain'}]; +% single interval up/down with cue +stimulusParameters.instructions{2}= [{'count the tones you hear clearly'}, { }, { 'ignore indistinct tones'}]; +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/multithreshold 1.46/paradigms/paradigm_thr_TMC.m Tue May 31 09:13:07 2011 +0100 @@ -0,0 +1,36 @@ +function paradigm_thr_TMC(handles) +global stimulusParameters experiment betweenRuns + +paradigm_training(handles) % default + +betweenRuns.variableName1='targetFrequency'; +betweenRuns.variableList1=1000; +betweenRuns.variableList1=str2num(get(handles.edittargetFrequency,'string')); +betweenRuns.variableName2='targetDuration'; +betweenRuns.variableList2=0.016; +betweenRuns.randomizeSequence=1; % 'random sequence' + +% delay > masker > gap > target + + +stimulusParameters.targetType='tone'; +stimulusParameters.targetPhase='sin'; +stimulusParameters.targetFrequency=1000; +stimulusParameters.targetDuration=0.016; +stimulusParameters.targetLevel=stimulusParameters.WRVstartValues(1); + +stimulusParameters.rampDuration=0.004; + +experiment.stopCriteria2IFC=[75 3 5]; +experiment.singleIntervalMaxTrials=[20]; + + +% forced choice window interval +stimulusParameters.AFCsilenceDuration=0.5; + +% instructions to user +% single interval up/down no cue +stimulusParameters.instructions{1}= [{'YES if you hear the tone clearly'}, { }, { 'NO if not (or you are uncertain'}]; +% single interval up/down with cue +stimulusParameters.instructions{2}= [{'count the tones you hear clearly'}, { }, { 'ignore indistinct tones'}]; +
--- a/multithreshold 1.46/paradigms/paradigm_training.m Fri May 27 17:50:40 2011 +0100 +++ b/multithreshold 1.46/paradigms/paradigm_training.m Tue May 31 09:13:07 2011 +0100 @@ -2,7 +2,8 @@ global stimulusParameters experiment betweenRuns stimulusParameters.subjectSampleRate=44100; % compatible with file input -stimulusParameters.subjectSampleRate=50000; % compatible with file input +% stimulusParameters.subjectSampleRate=64000; % compatible with file input +% stimulusParameters.subjectSampleRate=128000; % compatible with file input % assessment method % {'oneIntervalUpDown', 'MaxLikelihood', '2I2AFC++', '2I2AFC+++'}
--- a/multithreshold 1.46/printReport.m Fri May 27 17:50:40 2011 +0100 +++ b/multithreshold 1.46/printReport.m Tue May 31 09:13:07 2011 +0100 @@ -45,6 +45,7 @@ end end +clc fprintf('******** multithreshold ') x=pwd; disp(['version ' x(end-3:end)]) fprintf('\nName:\t%s ', experiment.name) @@ -191,10 +192,8 @@ fprintf('\n') end -% resultsSoFar=[betweenRuns.var1Sequence(betweenRuns.runNumber)'... -% betweenRuns.var2Sequence(betweenRuns.runNumber)'... -% betweenRuns.thresholds(betweenRuns.runNumber)' ]; -% fprintf('%10.3f \t%10.3f \t%10.1f \n', resultsSoFar') + +fprintf('\nparadigm:\t%s\n ', experiment.paradigm) % ------------------------------------------------------- sortTablesForPrinting function table= sortTablesForPrinting(idx1,idx2, var1values,var2values, x)
--- a/multithreshold 1.46/subjGUI_MT.m Fri May 27 17:50:40 2011 +0100 +++ b/multithreshold 1.46/subjGUI_MT.m Tue May 31 09:13:07 2011 +0100 @@ -1560,6 +1560,142 @@ end end +function [modelResponse, MacGregorResponse]=MAPmodel( MAPplot, method) + +global experiment stimulusParameters audio withinRuns +global outerMiddleEarParams DRNLParams AN_IHCsynapseParams + +savePath=path; +addpath('..\MAP') +modelResponse=[]; +MacGregorResponse=[]; + +% mono only (column vector) +audio=audio(:,1)'; + +% if stop button pressed earlier +if experiment.stop, return, end + +% -------------------------------------------------------------- run Model +MAPparamsName=experiment.name; +showPlotsAndDetails=experiment.MAPplot; +AN_spikesOrProbability='spikes'; + +% [response, method]=MAPsequenceSeg(audio, method, 1:8); +global ICoutput ANdt + MAP1_14(audio, 1/method.dt, method.nonlinCF,... + MAPparamsName, AN_spikesOrProbability); + +if showPlotsAndDetails + options.showModelParameters=0; + options.showModelOutput=1; + options.printFiringRates=1; + options.showACF=0; + options.showEfferent=1; + showMAP(options) +end + +% No response, probably caused by hitting 'stop' button +if isempty(ICoutput), return, end + +% MacGregor response is the sum total of all final stage spiking +MacGregorResponse= sum(ICoutput,1); % use IC + +% ---------------------------------------------------------- end model run + +dt=ANdt; +time=dt:dt:dt*length(MacGregorResponse); + +% group delay on unit response +MacGonsetDelay= 0.004; +MacGoffsetDelay= 0.022; + +% now find the response of the MacGregor model during the target presentation + group delay +switch experiment.threshEstMethod + case {'2I2AFC++', '2I2AFC+++'} + idx= time>stimulusParameters.testTargetBegins+MacGonsetDelay ... + & time<stimulusParameters.testTargetEnds+MacGoffsetDelay; + nSpikesTrueWindow=sum(MacGregorResponse(:,idx)); + idx=find(time>stimulusParameters.testNonTargetBegins+MacGonsetDelay ... + & time<stimulusParameters.testNonTargetEnds+MacGoffsetDelay); + nSpikesFalseWindow=sum(MacGregorResponse(:,idx)); + % nSpikesDuringTarget is +ve when more spikes are found + % in the target window + difference= nSpikesTrueWindow-nSpikesFalseWindow; + + if difference>0 + % hit + nSpikesDuringTarget=experiment.MacGThreshold+1; + elseif difference<0 + % miss (wrong choice) + nSpikesDuringTarget=experiment.MacGThreshold-1; + else + if rand>0.5 + % hit (random choice) + nSpikesDuringTarget=experiment.MacGThreshold+1; + else + % miss (random choice) + nSpikesDuringTarget=experiment.MacGThreshold-1; + end + end + disp(['level target dummy decision: ' ... + num2str([withinRuns.variableValue nSpikesTrueWindow ... + nSpikesFalseWindow nSpikesDuringTarget], '%4.0f') ] ) + + otherwise + % idx=find(time>stimulusParameters.testTargetBegins+MacGonsetDelay ... + % & time<stimulusParameters.testTargetEnds+MacGoffsetDelay); + % no delay at onset + idx=find(time>stimulusParameters.testTargetBegins +MacGonsetDelay... + & time<stimulusParameters.testTargetEnds+MacGoffsetDelay); + nSpikesDuringTarget=sum(MacGregorResponse(:,idx)); + + % find(MacGregorResponse)*dt-stimulusParameters.stimulusDelay + timeX=time(idx); +end + +% now find the response of the MacGregor model at the end of the masker +idx2=find(time>stimulusParameters.testTargetBegins-0.02 ... + & time<stimulusParameters.testTargetBegins); +if ~isempty(idx2) + maskerRate=mean(mean(MacGregorResponse(idx2))); +else + %e.g. no masker + maskerRate=0; +end + +if experiment.MAPplot + % add vertical lines to indicate target region + figure(99), subplot(6,1,6) + hold on + yL=get(gca,'YLim'); + plot([stimulusParameters.testTargetBegins + MacGonsetDelay ... + stimulusParameters.testTargetBegins + MacGonsetDelay],yL,'r') + plot([stimulusParameters.testTargetEnds + MacGoffsetDelay ... + stimulusParameters.testTargetEnds + MacGoffsetDelay],yL,'r') +end + +% specify unambiguous response +switch experiment.paradigm + case 'gapDetection' + gapResponse=(maskerRate-nSpikesDuringTarget)/maskerRate; + if gapResponse>0.2 + modelResponse=2; % gap detected + else + modelResponse=1; % gap not detected + end + [nSpikesDuringTarget maskerRate gapResponse modelResponse] + figure(22), plot(timeX,earObject(idx)) + otherwise + if nSpikesDuringTarget>experiment.MacGThreshold + modelResponse=2; % stimulus detected + else + modelResponse=1; % nothing heard (default) + end +end + + +path(savePath) % -----------------------------------------------------statsModelRunsGUI % The computer presses the buttons
--- a/multithreshold 1.46/testAN.m Fri May 27 17:50:40 2011 +0100 +++ b/multithreshold 1.46/testAN.m Tue May 31 09:13:07 2011 +0100 @@ -1,9 +1,9 @@ -function testAN +function vectorStrength=testAN(targetFrequency,BFlist, levels) % testIHC used either for IHC I/O function ... % or receptive field (doReceptiveFields=1) global experiment method stimulusParameters -global IHC_VResp_VivoParams IHCpreSynapseParams +global IHC_VResp_VivoParams IHC_cilia_RPParams IHCpreSynapseParams global AN_IHCsynapseParams % global saveMembranePotential MacGregorMultiParams dbstop if error @@ -12,8 +12,11 @@ ['..' filesep 'parameterStore'], ['..' filesep 'wavFileStore'],... ['..' filesep 'testPrograms']) -levels=-10:10:80; -% levels=80:10:90; +if nargin<3 + levels=-10:10:80; + % levels=80:10:90; +end + nLevels=length(levels); toneDuration=.2; @@ -22,8 +25,8 @@ localPSTHbinwidth=0.001; % Use only the first frequency in the GUI targetFrequency box to defineBF -targetFrequency=stimulusParameters.targetFrequency(1); -BFlist=targetFrequency; +% targetFrequency=stimulusParameters.targetFrequency(1); +% BFlist=targetFrequency; AN_HSRonset=zeros(nLevels,1); AN_HSRsaturated=zeros(nLevels,1); @@ -43,15 +46,23 @@ figure(15), clf set(gcf,'position',[980 356 401 321]) figure(5), clf - set(gcf,'position', [980 34 400 295]) +set(gcf,'position', [980 34 400 295]) +set(gcf,'name',[num2str(BFlist), ' Hz']); drawnow levelNo=0; for leveldB=levels levelNo=levelNo+1; - % sample rate should be amultiple of the targetFrequency for PSTH below + %% sample rate should be amultiple of the targetFrequency for PSTH below sampleRate=50000; + sampleRate=20*targetFrequency; + if sampleRate<20000 + sampleRate=round(40000/targetFrequency)*targetFrequency; + end + + + %% ananana dt=1/sampleRate; period=1/targetFrequency; dt=dt*(dt/period)*round(period/dt); @@ -79,8 +90,8 @@ global ANoutput CNoutput ICoutput ICmembraneOutput tauCas global ARattenuation MOCattenuation - - MAP1_14(inputSignal, 1/dt, targetFrequency, ... + + MAP1_14(inputSignal, 1/dt, BFlist, ... MAPparamsName, AN_spikesOrProbability); nTaus=length(tauCas); @@ -111,8 +122,8 @@ hold off, bar(PSTHtime,PSTH, 'b') hold on, bar(PSTHtime,PSTHLSR,'r') ylim([0 1000]) -xlim([0 length(PSTH)*localPSTHbinwidth]) - + xlim([0 length(PSTH)*localPSTHbinwidth]) + % AN - CV % CV is computed 5 times. Use the middle one (3) as most typical cvANHSR= UTIL_CV(HSRspikes, dt); @@ -180,7 +191,7 @@ xlim([0 max(time)]) title(['IC ' num2str(leveldB,'%4.0f') 'dB']) drawnow - + figure(5), subplot(2,2,1) plot(20*log10(MOC), 'k'), title(' MOC'), ylabel('dB attenuation') @@ -189,9 +200,9 @@ end % level figure(5), subplot(2,2,1) - plot(levels,20*log10(MOC), 'k'), - title(' MOC'), ylabel('dB attenuation') - ylim([-30 0]) +plot(levels,20*log10(MOC), 'k'), +title(' MOC'), ylabel('dB attenuation') +ylim([-30 0]) xlim([0 max(levels)]) fprintf('\n') @@ -209,7 +220,6 @@ % ---------------------------------------------------- display parameters figure(15), clf -set(gcf,'position',[1000 356 381 321]) nRows=2; nCols=2; % AN rate - level ONSET functions @@ -261,6 +271,7 @@ xlim([min(levels) max(levels)]) set(gca,'xtick',[levels(1):20:levels(end)]), grid on + ttl=['spont=' num2str(mean(ICHSRsaturated(1,:)),'%4.0f') ... ' sat=' num2str(mean(ICHSRsaturated(end,1)),'%4.0f')]; title( ttl) @@ -268,14 +279,13 @@ text(0, 350, 'IC', 'fontsize', 16) set(gcf,'name',' AN CN IC rate/level') -UTIL_showStruct(IHCpreSynapseParams, 'IHCpreSynapseParams') -UTIL_showStruct(AN_IHCsynapseParams, 'AN_IHCsynapseParams') +peakVectorStrength=max(vectorStrength); fprintf('\n') disp('levels vectorStrength') fprintf('%3.0f \t %6.4f \n', [levels; vectorStrength']) fprintf('\n') -fprintf('Phase locking, max vector strength= %6.4f\n\n',... +fprintf('Phase locking, max vector strength=\t %6.4f\n\n',... max(vectorStrength)) allData=[ levels' AN_HSRonset AN_HSRsaturated... @@ -287,3 +297,23 @@ fprintf('VS (phase locking)= \t%6.4f\n\n',... max(vectorStrength)) +UTIL_showStruct(IHC_cilia_RPParams, 'IHC_cilia_RPParams') +UTIL_showStruct(IHCpreSynapseParams, 'IHCpreSynapseParams') +UTIL_showStruct(AN_IHCsynapseParams, 'AN_IHCsynapseParams') + +fprintf('\n') +disp('levels vectorStrength') +fprintf('%3.0f \t %6.4f \n', [levels; vectorStrength']) +fprintf('\n') +fprintf('Phase locking, max vector strength= \t%6.4f\n\n',... + max(vectorStrength)) + +allData=[ levels' AN_HSRonset AN_HSRsaturated... + AN_LSRonset AN_LSRsaturated ... + CNHSRsaturated CNLSRrate... + ICHSRsaturated ICLSRsaturated]; +fprintf('\n levels \tANHSR Onset \tANHSR adapted\tANLSR Onset \tANLSR adapted\tCNHSR\tCNLSR\tICHSR \tICLSR \n'); +UTIL_printTabTable(round(allData)) +fprintf('VS (phase locking)= \t%6.4f\n\n',... + max(vectorStrength)) +
--- a/multithreshold 1.46/testBM.m Fri May 27 17:50:40 2011 +0100 +++ b/multithreshold 1.46/testBM.m Tue May 31 09:13:07 2011 +0100 @@ -19,7 +19,7 @@ % levels= 50; nLevels=length(levels); relativeFrequencies=[0.25 .5 .75 1 1.25 1.5 2]; -relativeFrequencies=1; +% relativeFrequencies=1; % refBMdisplacement is the displacement of the BM at threshold % 1 nm disp at threshold (9 kHz, Ruggero) @@ -157,4 +157,4 @@ UTIL_printTabTable([levels' finalSummary], ... num2str([0 stimulusFrequencies]','%5.0f'), '%5.0f') -path(savePath); \ No newline at end of file +path(savePath);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/multithreshold 1.46/testPhaseLocking.m Tue May 31 09:13:07 2011 +0100 @@ -0,0 +1,43 @@ +function testPhaseLocking + +testFrequencies=[250 500 1000 2000 4000 8000]; +levels=-20:10:80; +figure(14), clf +allStrengths=zeros(length(testFrequencies), length(levels)); +peakVectorStrength=zeros(1,length(testFrequencies)); +freqCount=0; +for targetFrequency=testFrequencies; + %single test + freqCount=freqCount+1; + vectorStrength=testAN(targetFrequency,targetFrequency, levels); + allStrengths(freqCount,:)=vectorStrength'; + peakVectorStrength(freqCount)=max(vectorStrength'); +end +%% plot results +figure(14) +subplot(2,1,2) +plot(levels,allStrengths) +xlabel('levels') +ylabel('vector strength') +legend (num2str(testFrequencies'),'location','eastOutside') + +subplot(2,1,1) +semilogx(testFrequencies,peakVectorStrength) +grid on +title ('peak vector strength') +xlabel('frequency') +ylabel('vector strength') + +johnson=[250 0.79 +500 0.82 +1000 0.8 +2000 0.7 +4000 0.25 +5500 0.05 +]; +hold on +plot(johnson(:,1),johnson(:,2),'o') +legend({'model','Johnson 80'},'location','eastOutside') +hold off + +
--- a/multithreshold 1.46/testRP.m Fri May 27 17:50:40 2011 +0100 +++ b/multithreshold 1.46/testRP.m Tue May 31 09:13:07 2011 +0100 @@ -9,7 +9,7 @@ figure(4), clf, set (gcf, 'name', ['IHC']) -% set(gcf,'position',[613 354 360 322]) +set(gcf,'position',[613 354 360 322]) drawColors='rgbkmcy'; drawnow @@ -204,4 +204,4 @@ UTIL_printTabTable([levels' allIHC_RP_dc]) path(savePath); -disp(paramChanges) \ No newline at end of file +disp(paramChanges)
--- a/multithreshold 1.46/testSynapse.m Fri May 27 17:50:40 2011 +0100 +++ b/multithreshold 1.46/testSynapse.m Tue May 31 09:13:07 2011 +0100 @@ -82,7 +82,7 @@ grid on figure(88), [c,H]=contour(time, maskerLevels,qtMatrix); clabel(c, H); -set(gcf,'position',[ 276 31 264 243]) +set(gcf,'position',[ 276 31 328 246]) xlabel('time'), ylabel('maskerLevels') path(savePath);
--- a/parameterStore/MAPparamsNormal.m Fri May 27 17:50:40 2011 +0100 +++ b/parameterStore/MAPparamsNormal.m Tue May 31 09:13:07 2011 +0100 @@ -80,7 +80,7 @@ % DRNL nonlinear path DRNLParams.a=3e4; % nonlinear path gain (below compression threshold) -% DRNLParams.a=0; % DRNL.a=0 means no OHCs (no nonlinear path) +% DRNLParams.a=3e2; % DRNL.a=0 means no OHCs (no nonlinear path) DRNLParams.b=8e-6; % *compression threshold raised compression % DRNLParams.b=1; % b=1 means no compression @@ -135,19 +135,19 @@ % #5 IHC_RP IHC_cilia_RPParams.Cab= 4e-012; % IHC capacitance (F) +IHC_cilia_RPParams.Cab= 1e-012; % IHC capacitance (F) IHC_cilia_RPParams.Et= 0.100; % endocochlear potential (V) -% IHC_cilia_RPParams.Et= 0.070; % endocochlear potential (V) +IHC_cilia_RPParams.Et= 0.07; % endocochlear potential (V) IHC_cilia_RPParams.Gk= 2e-008; % 1e-8 potassium conductance (S) IHC_cilia_RPParams.Ek= -0.08; % -0.084 K equilibrium potential IHC_cilia_RPParams.Rpc= 0.04; % combined resistances - %% #5 IHCpreSynapse IHCpreSynapseParams=[]; IHCpreSynapseParams.GmaxCa= 14e-9;% maximum calcium conductance -IHCpreSynapseParams.GmaxCa= 13e-9;% maximum calcium conductance +IHCpreSynapseParams.GmaxCa= 12e-9;% maximum calcium conductance IHCpreSynapseParams.ECa= 0.066; % calcium equilibrium potential IHCpreSynapseParams.beta= 400; % determine Ca channel opening IHCpreSynapseParams.gamma= 100; % determine Ca channel opening
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/testPrograms/demoTwisterProbability.m Tue May 31 09:13:07 2011 +0100 @@ -0,0 +1,130 @@ +function demoTwisterProbability + +% MAPdemo runs the MATLAB auditory periphery model (MAP1_14) as far as +% the AN (probabilities) or IC (spikes) with graphical output + +% Things you might want to change; #1 - #5 + +%% #1 parameter file name +MAPparamsName='Normal'; + + +%% #2 probability (fast) or spikes (slow) representation +% AN_spikesOrProbability='spikes'; +% or +AN_spikesOrProbability='probability'; + + +%% #3 pure tone, harmonic sequence or speech file input +signalType= 'tones'; +duration=0.100; % seconds +duration=0.020; % seconds +sampleRate= 64000; +% toneFrequency= 250:250:8000; % harmonic sequence (Hz) +toneFrequency= 2000; % or a pure tone (Hz8 + +rampDuration=.005; % seconds + +% or +signalType= 'file'; +fileName='twister_44kHz'; +% fileName='new-da-44khz'; + + +%% #4 rms level +% signal details +leveldBSPL=70; % dB SPL + + +%% #5 number of channels in the model +% 21-channel model (log spacing) +numChannels=21; +lowestBF=250; highestBF= 8000; +BFlist=round(logspace(log10(lowestBF), log10(highestBF), numChannels)); + +% or specify your own channel BFs +% BFlist=toneFrequency; + + +%% #6 change model parameters +paramChanges=[]; + +% or +% Parameter changes can be used to change one or more model parameters +% *after* the MAPparams file has been read +% This example declares only one fiber type with a calcium clearance time +% constant of 80e-6 s (HSR fiber) when the probability option is selected. +% switch AN_spikesOrProbability +% case 'probability' +% paramChanges={'IHCpreSynapseParams.tauCa=80e-6;'}; +% otherwise +% paramChanges=[]; +% end + +%% delare showMap options +showMapOptions=[]; % use defaults + +% or (example: show everything including an smoothed SACF output + showMapOptions.showModelParameters=1; + showMapOptions.showModelOutput=1; + showMapOptions.printFiringRates=1; + showMapOptions.showACF=0; + showMapOptions.showEfferent=1; + +%% Generate stimuli + +dbstop if error +restorePath=path; +addpath (['..' filesep 'MAP'], ['..' filesep 'wavFileStore']) +switch signalType + case 'tones' + inputSignal=createMultiTone(sampleRate, toneFrequency, ... + leveldBSPL, duration, rampDuration); + + case 'file' + [inputSignal sampleRate]=wavread(fileName); + inputSignal(:,1); + targetRMS=20e-6*10^(leveldBSPL/20); + rms=(mean(inputSignal.^2))^0.5; + amp=targetRMS/rms; + inputSignal=inputSignal*amp; +end + + +%% run the model +tic + +MAP1_14(inputSignal, sampleRate, BFlist, ... + MAPparamsName, AN_spikesOrProbability, paramChanges); +toc + +% the model run is now complete. Now display the results +showMAP(showMapOptions) + +toc +path(restorePath) + + +function inputSignal=createMultiTone(sampleRate, toneFrequency, ... + leveldBSPL, duration, rampDuration) +% Create pure tone stimulus +dt=1/sampleRate; % seconds +time=dt: dt: duration; +inputSignal=sum(sin(2*pi*toneFrequency'*time), 1); +amp=10^(leveldBSPL/20)*28e-6; % converts to Pascals (peak) +inputSignal=amp*inputSignal; + +% apply ramps +% catch rampTime error +if rampDuration>0.5*duration, rampDuration=duration/2; end +rampTime=dt:dt:rampDuration; +ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ... + ones(1,length(time)-length(rampTime))]; +inputSignal=inputSignal.*ramp; +ramp=fliplr(ramp); +inputSignal=inputSignal.*ramp; + +% add 10 ms silence +silence= zeros(1,round(0.03/dt)); +% inputSignal= [silence inputSignal silence]; +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/testPrograms/demoTwisterSpikes.m Tue May 31 09:13:07 2011 +0100 @@ -0,0 +1,134 @@ +function demoTwisterSpikes + +% MAPdemo runs the MATLAB auditory periphery model (MAP1_14) as far as +% the AN (probabilities) or IC (spikes) with graphical output + +% Things you might want to change; #1 - #5 + +%% #1 parameter file name +MAPparamsName='Normal'; + + +%% #2 probability (fast) or spikes (slow) representation +AN_spikesOrProbability='spikes'; +% or +% AN_spikesOrProbability='probability'; + + +%% #3 pure tone, harmonic sequence or speech file input +signalType= 'tones'; +duration=0.100; % seconds +duration=0.020; % seconds +sampleRate= 64000; +% toneFrequency= 250:250:8000; % harmonic sequence (Hz) +toneFrequency= 2000; % or a pure tone (Hz8 + +rampDuration=.005; % seconds + +% or +signalType= 'file'; +fileName='twister_44kHz'; +% fileName='new-da-44khz'; + + +%% #4 rms level +% signal details +leveldBSPL=70; % dB SPL + + +%% #5 number of channels in the model +% 21-channel model (log spacing) +numChannels=21; +lowestBF=250; highestBF= 8000; +BFlist=round(logspace(log10(lowestBF), log10(highestBF), numChannels)); + +% or specify your own channel BFs +% BFlist=toneFrequency; + + +%% #6 change model parameters +paramChanges=[]; + +% or +% Parameter changes can be used to change one or more model parameters +% *after* the MAPparams file has been read +% This example declares only one fiber type with a calcium clearance time +% constant of 80e-6 s (HSR fiber) when the probability option is selected. +% switch AN_spikesOrProbability +% case 'probability' +% paramChanges={'IHCpreSynapseParams.tauCa=80e-6;'}; +% otherwise +% paramChanges=[]; +% end + +%% delare showMap options +showMapOptions=[]; % use defaults + +% or (example: show everything including an smoothed SACF output + showMapOptions.showModelParameters=1; + showMapOptions.showModelOutput=1; + showMapOptions.printFiringRates=1; + showMapOptions.showACF=0; + showMapOptions.showEfferent=1; + +%% Generate stimuli + +dbstop if error +restorePath=path; +addpath (['..' filesep 'MAP'], ['..' filesep 'wavFileStore']) +switch signalType + case 'tones' + inputSignal=createMultiTone(sampleRate, toneFrequency, ... + leveldBSPL, duration, rampDuration); + + case 'file' + [inputSignal sampleRate]=wavread(fileName); + inputSignal(:,1); + targetRMS=20e-6*10^(leveldBSPL/20); + rms=(mean(inputSignal.^2))^0.5; + amp=targetRMS/rms; + inputSignal=inputSignal*amp; +end + + +%% run the model +tic + +fprintf('\n') +disp(['Signal duration= ' num2str(length(inputSignal)/sampleRate)]) +disp([num2str(numChannels) ' channel model']) +disp('Computing ...') +MAP1_14(inputSignal, sampleRate, BFlist, ... + MAPparamsName, AN_spikesOrProbability, paramChanges); +toc + +% the model run is now complete. Now display the results +showMAP(showMapOptions) + +toc +path(restorePath) + + +function inputSignal=createMultiTone(sampleRate, toneFrequency, ... + leveldBSPL, duration, rampDuration) +% Create pure tone stimulus +dt=1/sampleRate; % seconds +time=dt: dt: duration; +inputSignal=sum(sin(2*pi*toneFrequency'*time), 1); +amp=10^(leveldBSPL/20)*28e-6; % converts to Pascals (peak) +inputSignal=amp*inputSignal; + +% apply ramps +% catch rampTime error +if rampDuration>0.5*duration, rampDuration=duration/2; end +rampTime=dt:dt:rampDuration; +ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ... + ones(1,length(time)-length(rampTime))]; +inputSignal=inputSignal.*ramp; +ramp=fliplr(ramp); +inputSignal=inputSignal.*ramp; + +% add 10 ms silence +silence= zeros(1,round(0.03/dt)); +% inputSignal= [silence inputSignal silence]; +
--- a/testPrograms/test_MAP1_14.m Fri May 27 17:50:40 2011 +0100 +++ b/testPrograms/test_MAP1_14.m Tue May 31 09:13:07 2011 +0100 @@ -12,15 +12,16 @@ %% #2 probability (fast) or spikes (slow) representation AN_spikesOrProbability='spikes'; % or -AN_spikesOrProbability='probability'; +% AN_spikesOrProbability='probability'; %% #3 pure tone, harmonic sequence or speech file input signalType= 'tones'; duration=0.100; % seconds -sampleRate= 100000; +duration=0.020; % seconds +sampleRate= 64000; % toneFrequency= 250:250:8000; % harmonic sequence (Hz) -toneFrequency= 1000; % or a pure tone (Hz8 +toneFrequency= 2000; % or a pure tone (Hz8 rampDuration=.005; % seconds @@ -42,7 +43,7 @@ BFlist=round(logspace(log10(lowestBF), log10(highestBF), numChannels)); % or specify your own channel BFs -% BFlist=1000; +% BFlist=toneFrequency; %% #6 change model parameters @@ -53,12 +54,12 @@ % *after* the MAPparams file has been read % This example declares only one fiber type with a calcium clearance time % constant of 80e-6 s (HSR fiber) when the probability option is selected. -switch AN_spikesOrProbability - case 'probability' - paramChanges={'IHCpreSynapseParams.tauCa=80e-6;'}; - otherwise - paramChanges=[]; -end +% switch AN_spikesOrProbability +% case 'probability' +% paramChanges={'IHCpreSynapseParams.tauCa=80e-6;'}; +% otherwise +% paramChanges=[]; +% end %% delare showMap options showMapOptions=[]; % use defaults @@ -67,7 +68,7 @@ showMapOptions.showModelParameters=1; showMapOptions.showModelOutput=1; showMapOptions.printFiringRates=1; - showMapOptions.showACF=1; + showMapOptions.showACF=0; showMapOptions.showEfferent=1; %% Generate stimuli @@ -93,6 +94,10 @@ %% run the model tic +fprintf('\n') +disp(['Signal duration= ' num2str(length(inputSignal)/sampleRate)]) +disp([num2str(numChannels) ' channel model']) +disp('Computing ...') MAP1_14(inputSignal, sampleRate, BFlist, ... MAPparamsName, AN_spikesOrProbability, paramChanges); toc @@ -100,6 +105,7 @@ % the model run is now complete. Now display the results showMAP(showMapOptions) +toc path(restorePath) @@ -124,5 +130,5 @@ % add 10 ms silence silence= zeros(1,round(0.03/dt)); -inputSignal= [silence inputSignal silence]; +% inputSignal= [silence inputSignal silence];