diff userProgramsMathiasDietz/test_binaural.m @ 38:c2204b18f4a2 tip

End nov big change

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/userProgramsMathiasDietz/test_binaural.m	Mon Nov 28 13:34:28 2011 +0000
@@ -0,0 +1,433 @@
+function test_binaural
+% test_binaural is a first attempt to produce a binaural model
+%  incorporating MSO and IC models.
+% The monaural response is computed first for left and right stimuli
+%  before using the CN response as input to the binaural MSO model
+%  that, in turn, feeds a single cell IC model.
+%
+% The function has no arguments and everything is set up internally
+%
+% #1
+% Identify the file (in 'MAPparamsName') containing the model parameters
+%  the default is 'PL' which uses primary like neurons in the CN to
+%  simulate spherical bushy cells
+%
+% #2
+% Set AN_spikesOrProbability'). to 'spikes'
+%
+% #3
+% Choose between a tone signal or file input (in 'signalType')
+%
+% #4
+% Set the signal rms level (in leveldBSPL)
+%
+% #5
+% Identify the channels in terms of their best frequencies in the vector
+%  BFlist. This is currently a single-channel model, so only one BF needed
+%
+% #6
+% Last minute changes to the parameters fetched earlier can be made using
+%  the cell array of strings 'paramChanges'.
+%  Each string must have the same format as the corresponding line in the
+%  file identified in 'MAPparamsName'
+% Currently this is used to specify that only HSR fibers are used and
+%  for changing the current per AN spike at the CN dendrite
+%
+% #7
+% specify the parameters of the MSO cells in the MSOParams structure
+%
+% #8
+% specify the parameters of the IC cells in the ICMSOParams structure
+%
+% #9
+% identify the plots required from MAP1_14 (i.e. before the bonaural model)
+%
+% #10
+% Specify ITDs. The program cycles through different ITDs
+%
+
+global CNoutput dtSpikes
+dbstop if error
+restorePath=path;
+addpath (['..' filesep 'MAP'],    ['..' filesep 'wavFileStore'], ...
+    ['..' filesep 'utilities'])
+
+%%  #1 monaural model parameter file name
+MAPparamsName='PL';
+
+
+%% #2 'spikes' are mandatory for this model
+AN_spikesOrProbability='spikes';
+
+
+%% #3 pure tone, harmonic sequence or speech file input
+signalType= 'tones';
+sampleRate= 50000;
+duration=0.050;                 % seconds
+rampDuration=.005;              % raised cosine ramp (seconds)
+beginSilence=0.050;
+endSilence=0.050;
+toneFrequency= 750;            % or a pure tone (Hz)
+
+%   or
+% harmonic sequence (Hz)
+% F0=210;
+% toneFrequency= F0:F0:8000;
+
+%   or
+% signalType= 'file';
+% fileName='twister_44kHz';
+
+if strcmp(signalType, 'file')
+    % needed for labeling plot
+    showMapOptions.fileName=fileName;
+else
+    showMapOptions.fileName=[];
+end
+
+%% #4 rms level
+leveldBSPL= 70;
+
+
+%% #5 number of channels in the model
+BFlist=toneFrequency;
+
+%% #6 change model parameters
+
+paramChanges={'IHCpreSynapseParams.tauCa=80e-6;',...
+    'MacGregorMultiParams.currentPerSpike=0.800e-6;'};
+
+% Parameter changes can be used to change one or more model parameters
+%  *after* the MAPparams file has been read
+
+%% #7 MSO parameters
+MSOParams.nNeuronsPerBF=	10;   % N neurons per BF
+MSOParams.type = 'primary-like cell';
+MSOParams.fibersPerNeuron=4;   % N input fibers
+MSOParams.dendriteLPfreq=2000;  % dendritic filter
+MSOParams.currentPerSpike=0.11e-6; % (A) per spike
+MSOParams.currentPerSpike=0.5e-6; % (A) per spike
+MSOParams.Cap=4.55e-9;   % cell capacitance (Siemens)
+MSOParams.tauM=5e-4;     % membrane time constant (s)
+MSOParams.Ek=-0.01;      % K+ eq. potential (V)
+MSOParams.dGkSpike=3.64e-5; % K+ cond.shift on spike,S
+MSOParams.tauGk=	0.0012; % K+ conductance tau (s)
+MSOParams.Th0=	0.01;   % equilibrium threshold (V)
+MSOParams.c=	0.01;       % threshold shift on spike, (V)
+MSOParams.tauTh=	0.015;  % variable threshold tau
+MSOParams.Er=-0.06;      % resting potential (V)
+MSOParams.Eb=0.06;       % spike height (V)
+MSOParams.debugging=0;        % (special)
+MSOParams.wideband=0;         % special for wideband units
+
+%% #8 IC parameters
+ICchopperParams.nNeuronsPerBF=	10;   % N neurons per BF
+ICchopperParams.type = 'chopper cell';
+ICchopperParams.fibersPerNeuron=10;  % N input fibers
+ICchopperParams.dendriteLPfreq=50;   % dendritic filter
+ICchopperParams.currentPerSpike=50e-9; % *per spike
+ICchopperParams.currentPerSpike=100e-9; % *per spike
+ICchopperParams.Cap=1.67e-8; % ??cell capacitance (Siemens)
+ICchopperParams.tauM=0.002;  % membrane time constant (s)
+ICchopperParams.Ek=-0.01;    % K+ eq. potential (V)
+ICchopperParams.dGkSpike=1.33e-4; % K+ cond.shift on spike,S
+ICchopperParams.tauGk=	0.0005;% K+ conductance tau (s)
+ICchopperParams.Th0=	0.01; % equilibrium threshold (V)
+ICchopperParams.c=	0;        % threshold shift on spike, (V)
+ICchopperParams.tauTh=	0.02; % variable threshold tau
+ICchopperParams.Er=-0.06;    % resting potential (V)
+ICchopperParams.Eb=0.06;     % spike height (V)
+ICchopperParams.PSTHbinWidth=	1e-4;
+
+%% #9 delare 'showMap' options to control graphical output
+% this applies to the monaural input model only
+showMapOptions.printModelParameters=0;   % prints all parameters
+showMapOptions.showModelOutput=1;   % plot all stages if monaural input
+
+%% #10 ITDs
+% the program cycles through a range of stimulus ITDs
+ITDs=0e-6:100e-6:2000e-6;
+% ITDs=0; % single shot
+
+%% Now start computing!
+figure(98), clf, set(gcf, 'name', 'binaural demo')
+MSOcounts=zeros(1,length(ITDs));
+ICcounts=zeros(1,length(ITDs));
+ITDcount=0;
+for ITD=ITDs
+    ITDcount=ITDcount+1;
+    delaySamples=round(ITD* sampleRate);
+    %% Generate stimuli
+    switch signalType
+        case 'tones'
+            % 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 silence
+            intialSilence= zeros(1,round(beginSilence/dt));
+            finalSilence= zeros(1,round(endSilence/dt));
+            inputSignal= [intialSilence inputSignal finalSilence];
+
+        case 'file'
+            %% file input simple or mixed
+            [inputSignal sampleRate]=wavread(fileName);
+            dt=1/sampleRate;
+            inputSignal=inputSignal(:,1);
+            targetRMS=20e-6*10^(leveldBSPL/20);
+            rms=(mean(inputSignal.^2))^0.5;
+            amp=targetRMS/rms;
+            inputSignal=inputSignal*amp;
+            intialSilence= zeros(1,round(0.1/dt));
+            finalSilence= zeros(1,round(0.2/dt));
+            inputSignal= [intialSilence inputSignal' finalSilence];
+    end
+
+    %% run the monaural model twice
+    t=dt:dt:dt*length(inputSignal);
+    for ear={'left','right'}
+        figure(98), subplot(4,1,1), colour='b'; hold off
+        switch ear{1}
+            case 'right'
+                inputSignal=circshift(inputSignal', delaySamples)';
+                hold on
+                colour='r';
+        end
+        plot(t, inputSignal, colour)
+        title ('binaural inputs signals')
+        ylabel('Pa'), xlabel('time')
+        xlim([0 max(t)])
+
+        % call to monaural model
+        MAP1_14(inputSignal, sampleRate, BFlist, ...
+            MAPparamsName, AN_spikesOrProbability, paramChanges);
+
+        % the model run is now complete. Now display the results
+        UTIL_showMAP(showMapOptions, paramChanges)
+
+        % copy the CN inputspiking pattern to the binaural display
+        figure(98)
+        switch ear{1}
+            case 'left'
+                CNoutputLeft=CNoutput;
+                subplot(4,2,3)
+            case 'right'
+                CNoutputRight=CNoutput;
+                subplot(4,2,4)
+        end
+        plotInstructions=[];
+        plotInstructions.axes=gca;
+        plotInstructions.displaydt=dtSpikes;
+        plotInstructions.title= ['CN spikes ' ear{1}];
+        plotInstructions.rasterDotSize=2;
+        if sum(sum(CNoutput))<100
+            plotInstructions.rasterDotSize=3;
+        end
+        UTIL_plotMatrix(CNoutput, plotInstructions);
+
+    end % left/ right ear
+
+    %%  MSO
+    % run MSO model using left and right CN spikes
+    MSOspikes=MSO(CNoutputLeft,CNoutputRight, dtSpikes, MSOParams);
+    
+    sumspikes=sum(sum(MSOspikes));
+    disp(['ITD/ MSO spikes count= ' num2str([ITD sumspikes])])
+    MSOcounts(ITDcount)=sumspikes;
+    figure(98), subplot(4,2, 8), cla, hold off, plot(ITDs*1e6, MSOcounts)
+
+    %% IC chopper
+    % run IC model using all MSO spikes as input to a single IC cell
+    ICMSOspikes=ICchopper(MSOspikes, dtSpikes, ICchopperParams);
+    
+    sumspikes=sum(sum(ICMSOspikes));
+    disp(['ITD/ ICMSO spikes count= ' num2str([ITD sumspikes])])
+    ICcounts(ITDcount)=sumspikes;
+    figure(98), subplot(4,2,8),hold on, plot(ITDs*1e6, ICcounts,'r')
+    xlabel('ITD'), ylabel(' spike count')
+    title('MSO (blue)/ IC (red) spike counts')
+    legend({'MSO','IC'})
+
+end % ITDs
+
+path(restorePath)
+
+
+function MSOspikes=MSO(CNoutputLeft,CNoutputRight, dtSpikes, MSOparams)
+%%
+[nMSOcells nEpochs]=size(CNoutputLeft);
+inputCurrent=zeros(nMSOcells, nEpochs);
+MSOmembranePotential=inputCurrent;
+
+MSO_tauM=MSOparams.tauM;
+MSO_tauGk=MSOparams.tauGk;
+MSO_tauTh=MSOparams.tauTh;
+MSO_cap=MSOparams.Cap;
+MSO_c=MSOparams.c;
+MSO_b=MSOparams.dGkSpike;
+MSO_Th0=MSOparams.Th0;
+MSO_Ek=MSOparams.Ek;
+MSO_Eb= MSOparams.Eb;
+MSO_Er=MSOparams.Er;
+
+MSO_E=zeros(nMSOcells,1);
+MSO_Gk=zeros(nMSOcells,1);
+MSO_Th=MSO_Th0*ones(nMSOcells,1);
+
+% Dendritic filtering, all spikes are replaced by CNalphaFunction functions
+MSOdendriteLPfreq= MSOparams.dendriteLPfreq;
+MSOcurrentPerSpike=MSOparams.currentPerSpike;
+MSOspikeToCurrentTau=1/(2*pi*MSOdendriteLPfreq);
+t=dtSpikes:dtSpikes:5*MSOspikeToCurrentTau;
+MSO_CNalphaFunction= (MSOcurrentPerSpike / ...
+    MSOspikeToCurrentTau)*t.*exp(-t / MSOspikeToCurrentTau);
+% show alpha function
+% figure(84), subplot(4,2,2), plot(t,MSO_CNalphaFunction)
+% title(['LP cutoff ' num2str(MSOdendriteLPfreq)])
+
+% convert CN spikes to post-dendritic current
+CN_spikes=CNoutputLeft+CNoutputRight;
+for i=1:nMSOcells
+    x= conv2(CN_spikes(i,:),  MSO_CNalphaFunction);
+    inputCurrent(i,:)=x(1:nEpochs);
+end
+
+if MSO_c==0
+    % faster computation when threshold is stable (c==0)
+    for t=1:nEpochs
+        s=MSO_E>MSO_Th0;
+        dE = (-MSO_E/MSO_tauM + inputCurrent(:,t)/MSO_cap +...
+            (MSO_Gk/MSO_cap).*(MSO_Ek-MSO_E))*dtSpikes;
+        dGk=-MSO_Gk*dtSpikes/MSO_tauGk +MSO_b*s;
+        MSO_E=MSO_E+dE;
+        MSO_Gk=MSO_Gk+dGk;
+        MSOmembranePotential(:,t)=MSO_E+s.*(MSO_Eb-MSO_E)+MSO_Er;
+    end
+else
+    %  threshold is changing (MSO_c>0; e.g. bushy cell)
+    for t=1:nEpochs
+        dE = (-MSO_E/MSO_tauM + ...
+            inputCurrent(:,t)/MSO_cap + (MSO_Gk/MSO_cap)...
+            .*(MSO_Ek-MSO_E))*dtSpikes;
+        MSO_E=MSO_E+dE;
+        s=MSO_E>MSO_Th;
+        MSOmembranePotential(:,t)=MSO_E+s.*(MSO_Eb-MSO_E)+MSO_Er;
+        dGk=-MSO_Gk*dtSpikes/MSO_tauGk +MSO_b*s;
+        MSO_Gk=MSO_Gk+dGk;
+
+        % After a spike, the threshold is raised
+        % otherwise it settles to its baseline
+        dTh=-(MSO_Th-MSO_Th0)*dtSpikes/MSO_tauTh +s*MSO_c;
+        MSO_Th=MSO_Th+dTh;
+    end
+end
+
+figure(98),subplot(4,1,3)
+hold off, imagesc(MSOmembranePotential)
+title ('MSO (V)')
+
+MSOspikes=MSOmembranePotential> -0.01;
+% Remove any spike that is immediately followed by a spike
+%  NB 'find' works on columns (whence the transposing)
+MSOspikes=MSOspikes';
+idx=find(MSOspikes);
+idx=idx(1:end-1);
+MSOspikes(idx+1)=0;
+MSOspikes=MSOspikes';
+
+
+function ICMSOspikes=ICchopper(ICMSOspikes, dtSpikes, ICMSOParams)
+%%
+ICMSOspikes=sum(ICMSOspikes);
+[nICMSOcells nEpochs]=size(ICMSOspikes);
+inputCurrent=zeros(nICMSOcells, nEpochs);
+ICMSOmembranePotential=inputCurrent;
+
+ICMSO_tauM=ICMSOParams.tauM;
+ICMSO_tauGk=ICMSOParams.tauGk;
+ICMSO_tauTh=ICMSOParams.tauTh;
+ICMSO_cap=ICMSOParams.Cap;
+ICMSO_c=ICMSOParams.c;
+ICMSO_b=ICMSOParams.dGkSpike;
+ICMSO_Th0=ICMSOParams.Th0;
+ICMSO_Ek=ICMSOParams.Ek;
+ICMSO_Eb= ICMSOParams.Eb;
+ICMSO_Er=ICMSOParams.Er;
+
+ICMSO_E=zeros(nICMSOcells,1);
+ICMSO_Gk=zeros(nICMSOcells,1);
+ICMSO_Th=ICMSO_Th0*ones(nICMSOcells,1);
+
+% Dendritic filtering, all spikes are replaced by CNalphaFunction functions
+ICMSOdendriteLPfreq= ICMSOParams.dendriteLPfreq;
+ICMSOcurrentPerSpike=ICMSOParams.currentPerSpike;
+ICMSOspikeToCurrentTau=1/(2*pi*ICMSOdendriteLPfreq);
+t=dtSpikes:dtSpikes:5*ICMSOspikeToCurrentTau;
+ICMSOalphaFunction= (ICMSOcurrentPerSpike / ...
+    ICMSOspikeToCurrentTau)*t.*exp(-t / ICMSOspikeToCurrentTau);
+% show alpha function
+% figure(84), subplot(4,2,5), plot(t,ICMSOalphaFunction)
+% title(['IC MSO LP cutoff ' num2str(ICMSOdendriteLPfreq)])
+
+% post-dendritic current
+for i=1:nICMSOcells
+    x= conv2(1*ICMSOspikes(i,:),  ICMSOalphaFunction);
+    inputCurrent(i,:)=x(1:nEpochs);
+end
+
+if ICMSO_c==0
+    % faster computation when threshold is stable (c==0)
+    for t=1:nEpochs
+        s=ICMSO_E>ICMSO_Th0;
+        dE = (-ICMSO_E/ICMSO_tauM + inputCurrent(:,t)/ICMSO_cap +...
+            (ICMSO_Gk/ICMSO_cap).*(ICMSO_Ek-ICMSO_E))*dtSpikes;
+        dGk=-ICMSO_Gk*dtSpikes/ICMSO_tauGk +ICMSO_b*s;
+        ICMSO_E=ICMSO_E+dE;
+        ICMSO_Gk=ICMSO_Gk+dGk;
+        ICMSOmembranePotential(:,t)=ICMSO_E+s.*(ICMSO_Eb-ICMSO_E)+ICMSO_Er;
+    end
+else
+    %  threshold is changing (ICMSO_c>0; e.g. bushy cell)
+    for t=1:nEpochs
+        dE = (-ICMSO_E/ICMSO_tauM + ...
+            inputCurrent(:,t)/ICMSO_cap + (ICMSO_Gk/ICMSO_cap)...
+            .*(ICMSO_Ek-ICMSO_E))*dtSpikes;
+        ICMSO_E=ICMSO_E+dE;
+        s=ICMSO_E>ICMSO_Th;
+        ICMSOmembranePotential(:,t)=ICMSO_E+s.*(ICMSO_Eb-ICMSO_E)+ICMSO_Er;
+        dGk=-ICMSO_Gk*dtSpikes/ICMSO_tauGk +ICMSO_b*s;
+        ICMSO_Gk=ICMSO_Gk+dGk;
+
+        % After a spike, the threshold is raised
+        % otherwise it settles to its baseline
+        dTh=-(ICMSO_Th-ICMSO_Th0)*dtSpikes/ICMSO_tauTh +s*ICMSO_c;
+        ICMSO_Th=ICMSO_Th+dTh;
+    end
+end
+
+t=dtSpikes:dtSpikes:dtSpikes*length(ICMSOmembranePotential);
+figure(98),subplot(4,2,7)
+plot(t, ICMSOmembranePotential)
+ylim([-0.07 0]), xlim([0 max(t)])
+title('IC (V)')
+
+ICMSOspikes=ICMSOmembranePotential> -0.01;
+% now remove any spike that is immediately followed by a spike
+% NB 'find' works on columns (whence the transposing)
+ICMSOspikes=ICMSOspikes';
+idx=find(ICMSOspikes);
+idx=idx(1:end-1);
+ICMSOspikes(idx+1)=0;
+ICMSOspikes=ICMSOspikes';
+