Mercurial > hg > map
diff testPrograms/test2toneSuppression.m @ 38:c2204b18f4a2 tip
End nov big change
| author | Ray Meddis <rmeddis@essex.ac.uk> |
|---|---|
| date | Mon, 28 Nov 2011 13:34:28 +0000 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/testPrograms/test2toneSuppression.m Mon Nov 28 13:34:28 2011 +0000 @@ -0,0 +1,216 @@ +% 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;