view testPrograms/test2toneSuppression.m @ 38:c2204b18f4a2 tip

End nov big change

% function test2toneSuppression
%
% Demonstration of two-tone suppression
%
% A probe tone is played at a fixed level and a second tone is introduced
%  half way through the presentation. The response to the combined signal
%  is recorded and analysed.
%
% The second tone called the seeep tone is presented at a reange of
%  frequencies and levels. In a linear system we might expect the addition
%  of a second tone to increase the magnitude of the output. 
%  In fact, it often results in a reduction in the response. 
%  This is called two-tone suppression.
%
% The effect of the sweep tone is represented in a countour plot showing
%  both reductions and increases in response.
% The background colour in this plot is the response to the 
%  fixed tone alone


% Ruggero et al 1992 (Fig 2)
fixedToneFrequency=8600;
fixedToneLevelsdB=[45 :5: 90];
% fixedToneLevelsdB=[30];
fixedToneDuration=.050;		      % seconds
sweeptoneLevelsdB=[ 0 : 10:90];
nSweepToneFrequencies=21;
sampleSweepFrequency=10600;

global  ANprobRateOutput DRNLoutput
dbstop if error
restorePath=path;
addpath (['..' filesep 'MAP'],    ['..' filesep 'wavFileStore'], ...
    ['..' filesep 'utilities'])
figure(5), clf
figure(87), clf

nFixedToneLevels=length(fixedToneLevelsdB);

lowestSweepFrequency=fixedToneFrequency/6;
highestSweepFrequency=fixedToneFrequency*3;
sweepToneFrequencies=round(logspace(log10(lowestSweepFrequency), ...
    log10(highestSweepFrequency), nSweepToneFrequencies));

% key channels for snapshots
[a BFchannel]=min((sweepToneFrequencies-fixedToneFrequency).^2);
[a sampleChannel]=min((sweepToneFrequencies-sampleSweepFrequency).^2);

sampleRate= max(44100, 4*highestSweepFrequency);
dt=1/sampleRate; % seconds

startSilenceDuration=0.010;
startSilence= zeros(1,startSilenceDuration*sampleRate);

sweepStartPTR=...
    round((startSilenceDuration + fixedToneDuration/2)*sampleRate);

BF_BMresponse=zeros(length(sweepToneFrequencies), ...
    length(fixedToneLevelsdB), length(sweeptoneLevelsdB));

fixedTonedBCount=0;
for fixedTonedB=fixedToneLevelsdB
    fixedTonedBCount=fixedTonedBCount+1;

    BMpeakResponse= zeros(length(sweeptoneLevelsdB),length(sweepToneFrequencies));
    ANpeakResponse= zeros(length(sweeptoneLevelsdB),length(sweepToneFrequencies));
    sweepToneLevelCount=0;
    for sweepToneDB=sweeptoneLevelsdB
        sweepToneLevelCount=sweepToneLevelCount+1;
        suppFreqCount=0;
        for sweepToneFrequency=sweepToneFrequencies
            suppFreqCount=suppFreqCount+1;

            % fixedTone tone
            time1=dt: dt: fixedToneDuration;
            amp=10^(fixedTonedB/20)*28e-6;
            inputSignal=amp*sin(2*pi*fixedToneFrequency*time1);
            rampDuration=.005; rampTime=dt:dt:rampDuration;
            ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ones(1,length(time1)-length(rampTime))];
            inputSignal=inputSignal.*ramp;
            inputSignal=inputSignal.*fliplr(ramp);
            nsignalPoints=length(inputSignal);
            sweepStart=round(nsignalPoints/2);
            nSweepPoints=nsignalPoints-sweepStart;

            % sweepTone
            time2= dt: dt: dt*nSweepPoints;
            % B: tone on
            amp=10^(sweepToneDB/20)*28e-6;
            inputSignal2=amp*sin(2*pi*sweepToneFrequency*time2);
            rampDuration=.005; rampTime=dt:dt:rampDuration;
            ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ones(1,length(time2)-length(rampTime))];
            inputSignal2=inputSignal2.*ramp;
            inputSignal2=inputSignal2.*fliplr(ramp);
            % add tone and sweepTone components
            inputSignal(sweepStart+1:end)= inputSignal(sweepStart+1:end)+inputSignal2;

            inputSignal=...
                [startSilence inputSignal ];

            %% run MAP
                MAPparamsName='Normal';
                AN_spikesOrProbability='probability';
                % only use HSR fibers (NB no acoustic reflex)
                paramChanges={'IHCpreSynapseParams.tauCa=80e-6;'};
                paramChanges={'IHCpreSynapseParams.tauCa=80e-6;',...
                                'DRNLParams.g=0;'};
                MAP1_14(inputSignal, sampleRate, fixedToneFrequency, ...
                    MAPparamsName, AN_spikesOrProbability, paramChanges);

                % find toneAlone response
                if sweepToneLevelCount==1 && suppFreqCount==1
                    BMfixedToneAloneRate=...
                        mean(abs(DRNLoutput(sweepStartPTR:end)));
                    ANfixedToneAloneRate=...
                        mean(abs(ANprobRateOutput(sweepStartPTR:end)));
                end
                
                   BF_BMresponse(suppFreqCount,fixedTonedBCount, ...
                       sweepToneLevelCount)=...
                        mean(abs(DRNLoutput(sweepStartPTR:end)));

                BMpeakResponse(sweepToneLevelCount,suppFreqCount)=...
                    mean(abs(DRNLoutput(sweepStartPTR:end)))...
                    /BMfixedToneAloneRate;
                ANpeakResponse(sweepToneLevelCount,suppFreqCount)=...
                    mean(abs(ANprobRateOutput(sweepStartPTR:end)))...
                    /ANfixedToneAloneRate;
                disp(['F2: ', num2str([sweepToneFrequency sweepToneDB ...
                    BMpeakResponse(sweepToneLevelCount,suppFreqCount)...
                    ANpeakResponse(sweepToneLevelCount,suppFreqCount)])...
                    ' dB'])

                figure(5)
                time=dt:dt:dt*length(inputSignal);
                subplot(3,1,1), plot(time, inputSignal)
                title(['stimulus: Suppressor=' ...
                    num2str([sweepToneFrequency, sweepToneDB]) ' Hz/ dB'])

                time=dt:dt:dt*length(DRNLoutput);
                subplot(3,1,2), plot(time, DRNLoutput)
                xlim([0 fixedToneDuration])
                        ylim([0 inf])

                time=dt:dt:dt*length(ANprobRateOutput);
                subplot(3,1,2), plot(time, ANprobRateOutput)
                xlim([0 fixedToneDuration])
                        ylim([0 500])
                title(['ANresponse: fixedTone' num2str([fixedToneFrequency, fixedTonedB]) ' Hz/ dB'])

                subplot(3,2,5)
                contourf(sweepToneFrequencies,sweeptoneLevelsdB,BMpeakResponse, ...
                    [.1:.1:.9 1.05] )
                set(gca,'Xscale','log')
                set(gca,'Xtick', [1000  4000],'xticklabel',{'1000', '4000'})
                set(gcf, 'name',['fixedToneFrequency= ' num2str(fixedToneFrequency)])
                title(['BM' num2str(fixedTonedB) ' dB'])

                subplot(3,2,6)
                contourf(sweepToneFrequencies,sweeptoneLevelsdB,ANpeakResponse, ...
                    [.1:.1:.9 1.05] )
                set(gca,'Xscale','log')
                set(gca,'Xtick', [1000  4000],'xticklabel',{'1000', '4000'})
                set(gcf, 'name',['fixedToneFrequency= ' num2str(fixedToneFrequency)])
                title(['AN:' num2str(fixedTonedB) ' dB'])
                drawnow
        end
    end

    %% plot matrix

    figure (87)
    subplot(3, nFixedToneLevels, fixedTonedBCount)
    surf(sweepToneFrequencies,sweeptoneLevelsdB,BMpeakResponse)
    set(gca,'Xscale','log')
    zlabel('gain')
    xlim([lowestSweepFrequency highestSweepFrequency])
    ylim([min(sweeptoneLevelsdB) max(sweeptoneLevelsdB)])
    title('BM response')
    view([-11 52])

    subplot(3, nFixedToneLevels, nFixedToneLevels+fixedTonedBCount)
    contourf(sweepToneFrequencies,sweeptoneLevelsdB,BMpeakResponse, ...
        [.1:.5:.9 0.99 1.05] )
    hold on
    plot(fixedToneFrequency, fixedTonedB, 'or', 'markerfacecolor','w')
    set(gca,'Xscale','log')
    set(gca,'Xtick', [1000  5000],'xticklabel',{'1000', '5000'})
    ylabel('(BM) sweep level')
    xlabel('(BM) sweep freq')
    title(['fixed tone level=' num2str(fixedTonedB) ' dB'])
%     colorbar

    subplot(3, nFixedToneLevels, 2*nFixedToneLevels+fixedTonedBCount)
    contourf(sweepToneFrequencies,sweeptoneLevelsdB,ANpeakResponse, ...
        [.1:.5:.9 0.99 1.05] )
    hold on
    plot(fixedToneFrequency, fixedTonedB, 'or', 'markerfacecolor','w')
    set(gca,'Xscale','log')
    set(gca,'Xtick', [1000  5000],'xticklabel',{'1000', '5000'})
    ylabel('(AN) sweep level')
    xlabel('(AN) sweep freq')
    title(['fixed tone level=' num2str(fixedTonedB) ' dB'])
%     colorbar
    
end

% Ruggero fig 2 (probe tone level is x-axis, sweep tone level is y-axis
figure(1),semilogy(fixedToneLevelsdB,squeeze(BF_BMresponse(sampleChannel,:,:)))
ylim([-inf inf])
legend(num2str(sweeptoneLevelsdB'),'location','southeast')
xlabel('probe SPL')
ylabel ('displacement (m)')
title(['Probe ' num2str(fixedToneFrequency) ' Hz. Sweep ' ...
    num2str(sweepToneFrequencies(sampleChannel)) ' Hz.'])
path=restorePath;