annotate testPrograms/testRP.m @ 38:c2204b18f4a2 tip

End nov big change

function testRP(BFs,MAPparamsName,paramChanges)
% testIHC evaluates IHC I/O function
% multiple BFs can be used but only one is easier to interpret.
% e.g. testRP(1000,'Normal',{});

global experiment method inputStimulusParams
global stimulusParameters IHC_VResp_VivoParams IHC_cilia_RPParams
savePath=path;
addpath (['..' filesep 'utilities'],['..' filesep 'MAP'])
dbstop if error

figure(4), clf,
set (gcf, 'name', ['IHC'])
set(gcf,'position',[613   354   360   322])
drawColors='rgbkmcy';
drawnow

if nargin<3
    paramChanges=[]; 
end
if nargin<2
    MAPparamsName='Normal';
end
if nargin<3
    BFs=800; 
end

levels=-20:10:100;
nLevels=length(levels);
toneDuration=.05;
silenceDuration=.01;
sampleRate=50000;
dt=1/sampleRate;

allIHC_RP_peak=[];
allIHC_RP_dc=[];

for BFno=1:length(BFs)
    BF=BFs(BFno);
    targetFrequency=BF;
    % OR
    %Patuzzi and Sellick test (see ELP & AEM, 2006)
    % targetFrequency=100;

    IHC_RP_peak=zeros(nLevels,1);
    IHC_RP_min=zeros(nLevels,1);
    IHC_RP_dc=zeros(nLevels,1);

    time=dt:dt:toneDuration;

    rampDuration=0.004;
    rampTime=dt:dt:rampDuration;
    ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
        ones(1,length(time)-length(rampTime))];
    ramp=ramp.*fliplr(ramp);

    silence=zeros(1,round(silenceDuration/dt));

    toneStartptr=length(silence)+1;
    toneMidptr=toneStartptr+round(toneDuration/(2*dt)) -1;
    toneEndptr=toneStartptr+round(toneDuration/dt) -1;

    levelNo=0;
    for leveldB=levels
        levelNo=levelNo+1;
        % replicate at all levels
        amp=28e-6*10^(leveldB/20);

        %% create signal (leveldB/ frequency)
        inputSignal=amp*sin(2*pi*targetFrequency'*time);
        inputSignal= ramp.*inputSignal;
        inputSignal=[silence inputSignal silence];
        inputStimulusParams.sampleRate=1/dt;
%         global IHC_ciliaParams

        %% disable efferent for fast processing
        method.DRNLSave=1;
        method.IHC_cilia_RPSave=1;
        method.IHCpreSynapseSave=1;
        method.IHC_cilia_RPSave=1;
        method.segmentDuration=-1;
        moduleSequence=1:4;

        %% run the model
        global  DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV...
            IHCoutput
        AN_spikesOrProbability='probability';

        MAP1_14(inputSignal, sampleRate, BF, ...
            MAPparamsName, AN_spikesOrProbability, paramChanges);

        % DRNL
        DRNLoutput=DRNLoutput;
        DRNL_peak(levelNo,1)=max(DRNLoutput(toneMidptr:toneEndptr));
        DRNL_min(levelNo,1)=min(DRNLoutput(toneMidptr:toneEndptr));
        DRNL_dc(levelNo,1)=mean(DRNLoutput(toneMidptr:toneEndptr));

        % cilia
        IHC_ciliaData=IHC_cilia_output;
        IHC_ciliaData=IHC_ciliaData;
        IHC_cilia_peak(levelNo,1)=...
            max(IHC_ciliaData(toneMidptr:toneEndptr));
        IHC_cilia_min(levelNo,1)=...
            min(IHC_ciliaData(toneMidptr:toneEndptr));
        IHC_cilia_dc(levelNo,1)=...
            mean(IHC_ciliaData(toneMidptr:toneEndptr));

        % RP
        IHC_RPData=IHCoutput;
        IHC_RPData=IHC_RPData;
        IHC_RP_peak(levelNo,1)=...
            max(IHC_RPData(toneMidptr:toneEndptr));
        IHC_RP_min(levelNo,1)=...
            min(IHC_RPData(toneMidptr:toneEndptr));
        IHC_RP_dc(levelNo,1)=...
            mean(IHC_RPData(toneMidptr:toneEndptr));

    end % level


    disp(['parameter file was: ' MAPparamsName])
    fprintf('\n')

    %% plot DRNL
    subplot(2,2,1)
%     referenceDisp= 9e-9*1000/targetFrequency;
%     plot(levels,20*log10(DRNL_peak/referenceDisp), drawColors(BFno), ...
%         'linewidth',2), hold on
    referenceDisp=10e-9;
    plot(levels,20*log10(DRNL_peak/referenceDisp), drawColors(BFno), ...
        'linewidth',2), hold on
    title([' DRNL peak:  ' num2str(BFs) ' Hz'])
    ylabel ('log10DRNL(m)'), xlabel('dB SPL')
    xlim([min(levels) max(levels)]), ylim([-10 50])
    grid on

    %% plot cilia displacement
    figure(4)
    subplot(2,2,2)
    restingIHC_cilia=IHCrestingCiliaCond;
    plot(levels, IHC_cilia_peak,'k', 'linewidth',2), hold on
%     plot(levels, IHC_cilia_min,'r', 'linewidth',2)
%     hold on,
%     plot([min(levels) max(levels)], ...
%         [restingIHC_cilia restingIHC_cilia], 'g')
    title(' IHC apical cond.')
    ylabel ('IHCcilia(conductance)'), xlabel('dB SPL')
    xlim([min(levels) max(levels)])
    grid on

    %% plot receptor potentials
    figure(4)
    subplot(2,2,3)
    % RP I/O function min and max
    restingRP=IHC_RP_peak(1);
    toPlot= [fliplr(IHC_RP_min(:,1)') IHC_RP_peak(:,1)'];
    microPa=   28e-6*10.^(levels/20);
    microPa=[-fliplr(microPa) microPa];
    plot(microPa,toPlot, drawColors(BFno), 'linewidth',2)
    % ylim([0 300])

    %% Dallos and Harris data
    dallosx=[-0.9	-0.1	-0.001	0.001	0.01	0.9];
    dallosy=[-8	-7.8	-6.5	11	16.5	22]/1000 + restingRP;
    hold on, plot(dallosx,dallosy, 'o')
    plot([-1 1], [restingRP restingRP], 'r')
    title(' Dallos(86) data at 800 Hz')
    ylabel ('receptor potential(V)'), xlabel('Pa')
    ylim([-0.08 -0.02]), xlim([-1 1])
    grid on

    %% RP I/O function min and max
    figure(4)
    subplot(2,2,4)
    restingRP=IHC_RP_peak(1);
    peakRP=max(IHC_RP_peak);
    plot(levels, IHC_RP_peak,drawColors(BFno), 'linewidth',2)
    hold on
    plot(levels, IHC_RP_dc, [drawColors(BFno) ':'], 'linewidth',2)
    hold on,
    plot([min(levels) max(levels)], [restingRP restingRP], 'r')
    xlim([min(levels) max(levels)])
    % animal data
    sndLevel=[5	15	25	35	45	55	65	75];
RPanimal=restingRP+[0.5	2	4.6	5.8	6.4	7.2	8	10.2]/1000;
% could be misleading when restingRP changes
RPanimal=-0.060+[0.5	2	4.6	5.8	6.4	7.2	8	10.2]/1000;
hold on, plot(sndLevel,RPanimal,'o')

    grid on
    title(['Et= ' num2str(IHC_cilia_RPParams.Et) ':  RP data 7 kHz Patuzzi'])
    ylabel ('RP(V) peak and DC'), xlabel('dB SPL')
    ylim([-0.08 -0.04])
    allIHC_RP_peak=[allIHC_RP_peak IHC_RP_peak];
    allIHC_RP_dc=[allIHC_RP_dc IHC_RP_dc];

    fprintf('level\t peak\t DC\n')
    UTIL_printTabTable([levels' IHC_RP_peak IHC_RP_dc])
    % disp(['restingIHC_cilia= ' num2str(restingIHC_cilia)])
    fprintf('peakRP= \t%6.3f', peakRP)
    fprintf('\nrestingRP= \t%6.3f', restingRP)
    fprintf('\ndifference= \t%6.3f\n', (peakRP-restingRP))
    drawnow
end
% UTIL_showStruct(IHC_VResp_VivoParams, 'IHC_VResp_VivoParams')
UTIL_showStruct(IHC_cilia_RPParams, 'IHC_cilia_RPParams')
    fprintf('level\t peak\n')
    UTIL_printTabTable([levels' allIHC_RP_peak])
    fprintf('level\t DC\n')
    UTIL_printTabTable([levels' allIHC_RP_dc])

path(savePath);
disp(paramChanges)