THE AUDITORY MODELING TOOLBOX
Go to function
DATA_BAUMGARTNER2017looming - Results from Baumgartner et al. (2017)
Program code:
function varargout = data_baumgartner2017looming( varargin )
%DATA_BAUMGARTNER2017looming Results from Baumgartner et al. (2017)
%
% Usage: data = data_baumgartner2017looming(dataFlag,measure)
% data = data_baumgartner2017looming(fig)
%
% Output parameters:
% data : structure that contains either HRTFs (id, and Obj) or
% experimental results including raw and averaged data
% (rawData,data, opional meta information).
% Statisitcs (stat) and figure handles (fig) are
% provided if requested.
%
% DATA_BAUMGARTNER2017LOOMING provides individually measured HRTFs and
% experimental results from Baumgartner et al. (2017). Use the fig
% flag to obtain data shown in figures from Baumgartner et al. (2017).
%
% The dataFlag flag may be used to choose between HRTFs and various
% experimental results:
%
% 'hrtf' HRTFs used in all experiments.
% 'pretest' Behavioral results of pre-test.
% 'exp1' Behavioral or ERP results of Exp. I. Default.
% 'exp2' Behavioral results of Exp. II.
%
% The measure flag may be one of:
%
% 'judgment' Judgment of relative distance change (motion direction).
% Default.
% 'rt' Response time.
% 'erp' ERP magnitude measures.
%
% The fig flag may be one of:
%
% 'fig1b' Effect of spectral contrast manipulation according to
% factor C on magnitude responses of listener-specific
% stimuli of Exp. I (B, Top) as well as their frequency-specific
% and overall loudness changes relative to C = 1. Shaded
% areas denote ??1 standard error of the mean (SEM; N = 15).
% Note that changes in overall loudness oppose the intended
% effect of contrast switch.
% 'fig2' Behavioral responses were more consistent for sounds
% perceived as approaching compared to sounds perceived as
% receding if instantaneous spectral changes were presented
% in continuous stimuli. Mean behavioral responses in the
% 3?AFC motion discrimination task of Exp. I (fist figure; N = 15)
% and the 2?AFC motion discrimination task of Exp. II
% (second figure; N = 13). Results of Exp. II are separated
% between trials presenting instantaneous but continuous
% (Left; as in Exp. I) and discontinuous (Right) spectral
% contrast switches. Decreasing spectral contrast switches
% (orange lines) were predominantly perceived as approaching
% (orange triangles), increasing spectral contrast switches
% (green lines) as receding (green triangles), and constant
% spectral contrasts (no lines) as static (gray squares).
% Statistical analyses focused on these predominant response
% associations. Values reflect mean ??1 SEM.
% 'fig3a' Extracted N1 and P2 amplitudes evoked by stimulus onset.
% Error bars reflect SEM.
% 'fig3b' Extracted N1 and P2 amplitudes evoked by stimulus switch.
% Error bars reflect SEM.
% 'fig3c' Significant cluster in time and space that is distinctive
% between decreasing and increasing spectral contrasts. (no
% plotting)
%
% Additional flags may be:
%
% 'plot' Plot results as published.
% 'no_plot' No plots. Default.
% 'stat' Analyze and display inferential statistics.
% 'nostat' No statistics. Default.
% 'onset' To use onset ERPs.
% 'switch' To use switch-ERPs. Default.
%
%
% Requirements:
% -------------
%
% 1) SOFA API v0.4.3 or higher from http://sourceforge.net/projects/sofacoustics for Matlab (in e.g. thirdparty/SOFA)
%
% 2) Data in hrtf/baumgartner2017looming and auxdata/baumgartner2017looming (downloaded on the fly)
%
% 3) Statistics Toolbox for Matlab (for some fig)
%
% Examples:
% ---------
%
% To display results of Fig.1B :
%
% data_baumgartner2017looming('fig1b','plot');
%
% To display results of Fig.2 :
%
% data_baumgartner2017looming('fig2','plot');
%
% To display results of Fig.3A :
%
% data_baumgartner2017looming('fig3a','plot');
%
% To display results of Fig.3B :
%
% data_baumgartner2017looming('fig3b','plot');
%
% References:
% R. Baumgartner, D. K. Reed, B. Tóth, V. Best, P. Majdak, H. S. Colburn,
% and B. Shinn-Cunningham. Asymmetries in behavioral and neural responses
% to spectral cues demonstrate the generality of auditory looming bias.
% Proceedings of the National Academy of Sciences, 2017. [1]http ]
%
% References
%
% 1. http://www.pnas.org/content/early/2017/08/16/1703247114.abstract
%
%
% Url: http://amtoolbox.org/amt-1.3.0/doc/data/data_baumgartner2017looming.php
% #Author: Robert Baumgartner
% This file is licensed unter the GNU General Public License (GPL) either
% version 3 of the license, or any later version as published by the Free Software
% Foundation. Details of the GPLv3 can be found in the AMT directory "licences" and
% at <https://www.gnu.org/licenses/gpl-3.0.html>.
% You can redistribute this file and/or modify it under the terms of the GPLv3.
% This file is distributed without any warranty; without even the implied warranty
% of merchantability or fitness for a particular purpose.
% definput.import={'amt_cache'};
definput.flags.expirement = {'exp1','exp2','pretest','hrtf'};
definput.flags.measure = {'judgment','rt','erp'};
definput.flags.erp = {'switch','onset'};
definput.flags.plot = {'no_plot','plot'};
definput.flags.stat = {'nostat','stat'};
definput.flags.fig = {'nofig','fig1b','fig2','fig3a','fig3b','fig3c'};
[flags]=ltfatarghelper({},definput,lower(varargin));
if flags.do_nofig
%% HRTFs
if flags.do_hrtf
id = amt_load('baumgartner2017looming','subjects.mat');
out = struct('id',{},'Obj',{});
for ii = 1:length(id.subjects)
out(ii).id = id.subjects{ii};
out(ii).Obj = SOFAload(fullfile(SOFAdbPath,'baumgartner2017looming',...
[id.subjects{ii},'_eq.sofa']));
end
varargout{1} = out;
return
end
%% Behavioral Results
if flags.do_judgment || flags.do_rt
raw=amt_load('baumgartner2017looming',[flags.expirement,'.mat']);
Nsubj = length(raw.data);
ISI = unique(raw.data(1).ISI);
nISI = length(ISI);
ISIlbl = {'continuous','discont.'};
RTprctile = 25;
conditions = {[0,1];[0,0.5];[0.5,1];[1,0];[0.5,0];[1,0.5];[0,0];[0.5,0.5];[1,1]};
condLabelPlot = { '0\leftrightarrow1';'0\leftrightarrow.5';'.5\leftrightarrow1';...
'C_1=C_2'};
response = {'receding','approaching','constant'};
%% Evaluate Data
pResp = nan(Nsubj,length(conditions),length(response),nISI);
RTall = pResp;
for ss = 1:Nsubj
E = raw.data(ss).judgment;
C12 = raw.data(ss).C12;
sISI = raw.data(ss).ISI;
RT = raw.data(ss).RT;
for iISI = 1:nISI
for cc = 1:length(conditions)
idc = C12(:,1) == conditions{cc}(1) & C12(:,2) == conditions{cc}(2) & sISI(:) == ISI(iISI);
N = sum(idc);
Irecede = E(idc) > 0;
pResp(ss,cc,1,iISI) = sum(Irecede) / N;
RTall(ss,cc,1,iISI) = percentile(RT(Irecede),RTprctile);
Iappr = E(idc) < 0;
pResp(ss,cc,2,iISI) = sum(Iappr) / N;
RTall(ss,cc,2,iISI) = percentile(RT(Iappr),RTprctile);
Iconst = E(idc) == 0;
pResp(ss,cc,3,iISI) = sum(Iconst) / N;
RTall(ss,cc,3,iISI) = percentile(RT(Iconst),RTprctile);
end
end
end
if flags.do_judgment
meas = 100*pResp; % responses in percent
YLabel = 'Response (%)';
if flags.do_exp2
YLim = [43,105];
else
YLim = [-5,105];
end
elseif flags.do_rt
meas = 1e3*RTall; % response times in ms
YLabel = 'Response time (ms)';
YLim = 1e3*[0.451,1.149];
end
% outlier removal
rawData = struct2cell(raw.data);
ID = rawData(1,:);
if flags.do_exp2
pc = reshape(cat(2,pResp(:,1:3,1,:),pResp(:,4:6,2,:)),[Nsubj,6,nISI]); % percent correct
bias = squeeze(mean(pc(:,4:6,:) - pc(:,1:3,:),2));
[~,outlier] = max(bias(:,2)); % 1 listener showed a markedly larger looming bias for discontinuous stimuli compared with continuous stimuli (data provided in Outlier Evaluation for Experiment II and Fig. S3)
amt_disp(['Subject ' raw.data(outlier).ID ' identified as outlier and removed from further analyses.']);
iNew = (1:Nsubj)~=outlier;
ID = ID(iNew);
meas = meas(iNew,:,:,:);
Nsubj = Nsubj-1;
end
% Average constant conditions
meas = cat(2,meas(:,1:6,:,:),nan_mean(meas(:,7:9,:,:),2));
conditions = [conditions(1:6);'constant'];
% Standard errors of the means
sem = nan_std(meas)/sqrt(Nsubj);
% Output
out.data = meas;
out.rawData = raw.data;
out.meta.dim = 'subject_C_response_ISI';
out.meta.subject = ID;
out.meta.C = conditions;
out.meta.response = response;
out.meta.ISI = ISI;
if flags.do_stat && ~verLessThan('matlab','8.2')
% ANOVA
pc = reshape(cat(2,meas(:,1:3,1,:),meas(:,4:6,2,:)),Nsubj,[]); % percent correct
out.pcorrect.data = pc;
DV = array2table(pc);
contrast = strrep(condLabelPlot(1:3),'\leftrightarrow','-');
contrast = repmat(contrast,[2*nISI,1]);
out.pcorrect.contrast = contrast;
direction = cell(6,1);
direction(1:3) = {'receding'};
direction(4:6) = {'approaching'};
direction = repmat(direction,[nISI,1]);
out.pcorrect.direction = direction;
idISI = ceil((1:6*nISI)/6);
ISInom = nominal(ISI(idISI))';
IVs = table(contrast,direction,ISInom);
rm = fitrm(DV,['pc1-pc',num2str(length(contrast)),' ~ 1'],'WithinDesign',IVs);
if length(ISI) == 1
[ranovaResult,~,C,~] = ranova(rm,'WithinModel','direction*contrast');
else
[ranovaResult,~,C,~] = ranova(rm,'WithinModel','direction*contrast*ISInom');
end
ranovaResult.Properties.RowNames = strrep(ranovaResult.Properties.RowNames,'(Intercept):','');
% Sphericity corrections
spherCorr = epsilon(rm,C);
% Add corrected DFs to ranova table
idrep = round(0.5:0.5:length(spherCorr.GreenhouseGeisser)); % repeat iteratively
ranovaResult.DFGG = ranovaResult.DF .* ...
reshape(spherCorr.GreenhouseGeisser(idrep),size(ranovaResult.DF));
% Add effect sizes to ranova table
SSeffect = ranovaResult.SumSq(1:2:end);
SSerror = ranovaResult.SumSq(2:2:end);
eta_pSq = nan(2*length(SSerror),1);
eta_pSq(1:2:end) = SSeffect./(SSeffect+SSerror); % effect size per (eta_partial)^2
ranovaResult.eta_pSq = eta_pSq;
amt_disp(ranovaResult(:,[4,6,9,10]),'documentation');
mc = multcompare(rm,'contrast');
amt_disp(mc,'documentation');
out.stat = ranovaResult;
end
%% Plot
if flags.do_plot
out.fig = figure;
x = [1:3,1:3,4];
XTick = [x(1:3),x(end)];
if flags.do_exp2
x = x(1:6);
XTick = x(1:3);
response = response(1:2);
end
XLim = [XTick(1)-.6,XTick(end)+.6];
dx2 = .05*[-1,0,1]; % between responses
symb = '^vs';
lineStyle = {':';'-';'-.'}; % increase, decrease
colorR = [5,120,100;250,30,0;149,123,109]/255;
colorC = [5,120,100;250,30,0;255,255,255]/255;
for iisi = 1:nISI
subplot(1,nISI,iisi)
for rr = 1:length(response)
y = nan_mean(meas(:,:,rr,iisi));
l = sem(1,:,rr,iisi);
u = sem(1,:,rr,iisi);
idx = {1:3;4:6;7};
if flags.do_exp2
idx = idx(1:2);
end
for ii = 1:length(idx)
hC(rr,ii) = plot(x(idx{ii})+dx2(rr),y(idx{ii}),lineStyle{ii});
hold on
hR(rr,ii) = errorbar(x(idx{ii})+dx2(rr),y(idx{ii}),l(idx{ii}),u(idx{ii}),symb(rr));
set(hC(rr,ii),'Color',colorC(ii,:))
end
set(hR(rr,:),'MarkerFaceColor',colorR(rr,:),'Color',colorR(rr,:))
end
set(hR(1:2:size(hR,1),:),'MarkerFaceColor','w')
if flags.do_exp2
set([hC(2:3),hR(2:3)],'Visible','off')
end
set(gca,'XTick',XTick,'XTickLabel',condLabelPlot(1:length(XTick)))
axis([XLim,YLim])
ylabel(YLabel)
xlabel('Spectral contrast pair')
if flags.do_exp2
title(['Exp. II: ',ISIlbl{iisi}])
end
end
stimulus = {'C increase','C decrease','C constat'};
legend([hC(1,:)';hR(:,1)],[stimulus(1:size(hC,2)),response],'Location','eastoutside')
end
end
%% Physiological Results
if flags.do_erp
if not(flags.do_exp1)
error('ERPs only available for exp1.')
end
if flags.do_onset % onset
erp = amt_load('baumgartner2017looming','onsetERP.mat');
else % flags.do_switch
erp = amt_load('baumgartner2017looming','switchERP.mat');
end
erp.compLbl = {'N1','P2'};
out.rawData = erp;
if flags.do_stat
for rr = 1:length(erp.compLbl)
amt_disp(erp.compLbl{rr},'documentation');
amt_disp(erp.compStats{rr}.ranova(:,[4,6,9,10]),'documentation');
if isfield(erp.compStats{rr},'posthoc')
if isfield(erp.compStats{rr}.posthoc,'combination')
amt_disp(erp.compStats{rr}.posthoc.combination,'documentation');
else
amt_disp(erp.compStats{rr}.posthoc,'documentation');
end
end
end
end
if flags.do_plot
if flags.do_onset % onset
condLabelData = {'0','0.5','1'};
condLabel = condLabelData;
legLbl = erp.compLbl;
XLabel = 'Spectral contrast';
YLim = [-3.8,4.9];
condOrder = [1,3,2]; % to reorder erp.condLbl acc. to condLabel
resp = permute(erp.compAmp(condOrder,:,:),[2,1,3]); % subjects in first dimension
idx = {1:3};
dx = 0;
else % flags.do_switch
condLabel = { '0\leftrightarrow1';'0\leftrightarrow.5';'.5\leftrightarrow1'};
legLbl = {[erp.compLbl{1},', C increase'],[erp.compLbl{1},', C decrease'],...
[erp.compLbl{2},', C increase'],[erp.compLbl{2},', C decrease']};
XLabel = 'Spectral contrast pair';
YLim = [-3.5,3.5];
condOrder = [5,6,4,2,3,1]; % to reorder erp.condLbl acc. to condLabel
resp = permute(erp.compAmp(condOrder,:,:),[2,1,3]);
idx = {1:3;4:6};
dx = .1*[-1,1];
end
% Standard errors
seResp = std(resp)/sqrt(size(erp.compAmp,2));
out.fig = figure;
x = 1:3;
symb = 'oo';
lineStyle = {':','-'};
color = 0.8*[.65,.35,1;1,.5,0];
for rr = 1:length(erp.compLbl)
y = mean(resp(:,:,rr));
l = seResp(1,:,rr);
u = seResp(1,:,rr);
for ii = 1:length(idx)
h(rr,ii) = errorbar(x+dx(ii),y(idx{ii}),l(idx{ii}),u(idx{ii}),...
[symb(rr),lineStyle{ii}]);
hold on
end
set(h(rr,1),'MarkerFaceColor','w','Color',color(rr,:))
set(h(rr,length(idx)),'MarkerFaceColor',color(rr,:),'Color',color(rr,:))
set(gca,'XTick',x,'XLim',x([1,3])+[-1,1])
set(gca,'XTickLabel',condLabel)
ylabel('Cz potential (uV)')
xlabel(XLabel)
end
set(gca,'YLim',YLim,'YMinorTick','on')
legend(legLbl,'Location','eastoutside')
end
end
end
%% Fig. 1B
if flags.do_fig1b
stim = sig_baumgartner2017looming( 'exp1');
%% Top panel: Transfer characterisitcs
fs = stim(1).fs;
Nfft = 2^10;
freq = 0:fs/Nfft:fs/2; % frequency vector
mag = nan(length(freq),length(stim(1).C_IR),2,length(stim));
for ss = 1:length(stim)
if stim(ss).azi == -90
ipsiContra = [2,1];
else
ipsiContra = [1,2];
end
for cc = 1:length(stim(1).C_IR)
Sig = stim(ss).IR{cc};
for ich = 1:2
ch = ipsiContra(ich);
mag(:,cc,ich,ss) = db(abs(fftreal(Sig(:,ch),Nfft)));
end
end
end
mag = mag-3; % arbitrary adjustment to set ipsi. C=0 at 0 dB
if flags.do_plot
XLim = [800,17e3];
YLim = [-34,12];
blue = [0,0,0.7];
green = [0,0.7,0];
red = [0.7,0,0];
color = {blue,1.4*blue;green,1.4*green;red,1.4*red};
lineStyle = {'-','--'};
out.fig(1) = figure;
ii = 1;
for cc = 1:3
for ch = 1:2
lMEAN = mean(mag(:,cc,ch,:),4);
lSEM = std(mag(:,cc,ch,:),0,4)/sqrt(15);
h(ii) = shadedErrorBar(freq,lMEAN,lSEM,{'LineStyle',lineStyle{ch},'Color',color{cc,ch}},1);
hold on
ii = ii+1;
end
end
set(gca,'XScale','log','XLim',XLim,'YLim',YLim)
leg = legend([h.mainLine],'C=0 (ipsi)','C=0 (contra)','C=0.5 (ipsi)','C=0.5 (contra)','C=1 (ipsi)','C=1 (contra)');
set(leg,'Location','eastoutside','box','off')
ylabel('Magnitude (dB)')
xlabel('Frequency (Hz)')
end
out.magnitudeResponse.data = permute(mag,[1,3,4,2]);
out.magnitudeResponse.meta.dim = 'freq_channel_subject_C';
out.magnitudeResponse.meta.freq = freq;
out.magnitudeResponse.meta.channel = {'ipsi','contra'};
out.magnitudeResponse.meta.subject = cat(1,stim.ID);
out.magnitudeResponse.meta.C = 0:0.5:1;
%% Bottom panels: Loudness predictions
% The following loudness predictions were performed with the
% LoudnessToolbox 1.2 provided by Genesis (http://genesis-acoustics.com);
% models used:
% M1: Loudness_ISO532B_from_sound (calculated but not shown)
% M2: Loudness_ANSI_S34_2007 (used for publication);
% Data dimensions:
% subject (1:15) x C (0:.5:1) x model (M1,M2) [x channel (ipsi,contra)];
% frequency (M1:BarkAxis; M2:fc) x channel (ipsi,contra)
L = amt_load('baumgartner2017looming','specificLoudness.mat');
% Difference to reference
dLL_specif = cell(3,1);
for ss = 1:length(stim)
LLC1 = L.loudnessLevel_specif{ss,3,2};
for m = 1:3
dLL_specif{m}(:,:,ss) = L.loudnessLevel_specif{ss,m,2} - LLC1;
end
end
dLoudnessLevel = L.loudnessLevel(:,1:2,:,:) - repmat(L.loudnessLevel(:,3,:,:),[1,2,1,1]);
if flags.do_plot
out.fig(2) = figure;
YLim = [-9,13];
ii = 1;
for m = 1:3
for ch = 1:2
lMEAN = mean(dLL_specif{m}(:,ch,:),3);
lSEM = std(dLL_specif{m}(:,ch,:),0,3)/sqrt(15);
h(ii) = shadedErrorBar(L.fc,lMEAN,lSEM,{'LineStyle',lineStyle{ch},'Color',color{m,ch}},1);
hold on
ii = ii+1;
end
end
set(gca,'XScale','log','XLim',XLim,'YLim',YLim)
leg = legend([h.mainLine],'C=0 (ipsi)','C=0 (contra)','C=.5 (ipsi)','C=.5 (contra)');
set(leg,'Location','northwest','box','off')
ylabel('Loudness level difference to C=1 (phon)')
xlabel('Frequency (Hz)')
end
out.specificLoudnessLevelDiff.data = cat(4,dLL_specif{:});
out.specificLoudnessLevelDiff.meta.dim = 'freq_channel_subject_C';
out.specificLoudnessLevelDiff.meta.freq = L.fc;
out.specificLoudnessLevelDiff.meta.channel = {'ipsi','contra'};
out.specificLoudnessLevelDiff.meta.subject = cat(1,stim.ID);
out.specificLoudnessLevelDiff.meta.C = 0:0.5:1;
% overall loudness level
dLoudnessLevelP = dLoudnessLevel(:,:,2,:);
if flags.do_plot
try
out.fig(3) = figure;
boxplot(dLoudnessLevelP(:,:),... %,{{'M1';'M1';'M2';'M2'},[0;0.5;0;.5]}
'Factorgap',10,'FactorSeparator',1,'Whisker',Inf,...
'Colors',cat(1,color{[1,2,4,5]}))
set(gca,'YLim',YLim)
ylabel('Loudness level difference to C=1 (phon)')
xlabel('Spectral contrast (C)')
catch
warning('Statistics Toolbox not available, omitting figure 3.')
end
end
out.loudnessLevelDiff.data = permute(dLoudnessLevelP,[4,1,2,3]);
out.loudnessLevelDiff.meta.dim = 'channel_subject_C';
out.loudnessLevelDiff.meta.channel = {'ipsi','contra'};
out.loudnessLevelDiff.meta.subject = cat(1,stim.ID);
out.loudnessLevelDiff.meta.C = [0,0.5];
end
%% Fig. 2
if flags.do_fig2
amt_disp('Exp. I:','documentation');
out.exp1 = data_baumgartner2017looming('exp1',flags.plot,'stat');
title('Exp. I')
amt_disp('Exp. II:','documentation');
out.exp2 = data_baumgartner2017looming('exp2',flags.plot,'stat');
legend off
end
%% Fig. 3A
if flags.do_fig3a
out = data_baumgartner2017looming('erp','onset',flags.plot,flags.stat);
end
%% Fig. 3B
if flags.do_fig3b
out = data_baumgartner2017looming('erp','switch',flags.plot,flags.stat);
end
%% Fig. 3C
if flags.do_fig3c
out = amt_load('baumgartner2017looming','ERPclusterAnalysis.mat');
end
%% Output
if nargout == 1
varargout{1} = out;
end
end
%% Internal plotting functions
function varargout=shadedErrorBar(x,y,errBar,lineProps,transparent)
% function H=shadedErrorBar(x,y,errBar,lineProps,transparent)
%
% Purpose
% Makes a 2-d line plot with a pretty shaded error bar made
% using patch. Error bar color is chosen automatically.
%
% Inputs
% x - vector of x values [optional, can be left empty]
% y - vector of y values or a matrix of n observations by m cases
% where m has length(x);
% errBar - if a vector we draw symmetric errorbars. If it has a size
% of [2,length(x)] then we draw asymmetric error bars with
% row 1 being the upper bar and row 2 being the lower bar
% (with respect to y). ** alternatively ** errBar can be a
% cellArray of two function handles. The first defines which
% statistic the line should be and the second defines the
% error bar.
% lineProps - [optional,'-k' by default] defines the properties of
% the data line. e.g.:
% 'or-', or {'-or','markerfacecolor',[1,0.2,0.2]}
% transparent - [optional, 0 by default] if ==1 the shaded error
% bar is made transparent, which forces the renderer
% to be openGl. However, if this is saved as .eps the
% resulting file will contain a raster not a vector
% image.
%
% Outputs
% H - a structure of handles to the generated plot objects.
%
%
% Examples
% y=randn(30,80); x=1:size(y,2);
% shadedErrorBar(x,mean(y,1),std(y),'g');
% shadedErrorBar(x,y,{@median,@std},{'r-o','markerfacecolor','r'});
% shadedErrorBar([],y,{@median,@std},{'r-o','markerfacecolor','r'});
%
% Overlay two transparent lines
% y=randn(30,80)*10; x=(1:size(y,2))-40;
% shadedErrorBar(x,y,{@mean,@std},'-r',1);
% hold on
% y=ones(30,1)*x; y=y+0.06*y.^2+randn(size(y))*10;
% shadedErrorBar(x,y,{@mean,@std},'-b',1);
% hold off
%
%
% Rob Campbell - November 2009
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Error checking
narginchk(3,5)
%Process y using function handles if needed to make the error bar
%dynamically
if iscell(errBar)
fun1=errBar{1};
fun2=errBar{2};
errBar=fun2(y);
y=fun1(y);
else
y=y(:)';
end
if isempty(x)
x=1:length(y);
else
x=x(:)';
end
%Make upper and lower error bars if only one was specified
if length(errBar)==length(errBar(:))
errBar=repmat(errBar(:)',2,1);
else
s=size(errBar);
f=find(s==2);
if isempty(f), error('errBar has the wrong size'), end
if f==2, errBar=errBar'; end
end
if length(x) ~= length(errBar)
error('length(x) must equal length(errBar)')
end
%Set default options
defaultProps={'-k'};
if nargin<4, lineProps=defaultProps; end
if isempty(lineProps), lineProps=defaultProps; end
if ~iscell(lineProps), lineProps={lineProps}; end
if nargin<5, transparent=0; end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot to get the parameters of the line
H.mainLine=plot(x,y,lineProps{:});
% Work out the color of the shaded region and associated lines
% Using alpha requires the render to be openGL and so you can't
% save a vector image. On the other hand, you need alpha if you're
% overlaying lines. There we have the option of choosing alpha or a
% de-saturated solid colour for the patch surface .
col=get(H.mainLine,'color');
edgeColor=col+(1-col)*0.55;
patchSaturation=0.15; %How de-saturated or transparent to make patch
if transparent
faceAlpha=patchSaturation;
patchColor=col;
set(gcf,'renderer','openGL')
else
faceAlpha=1;
patchColor=col+(1-col)*(1-patchSaturation);
set(gcf,'renderer','painters')
end
%Calculate the error bars
uE=y+errBar(1,:);
lE=y-errBar(2,:);
%Add the patch error bar
holdStatus=ishold;
if ~holdStatus, hold on, end
%Make the patch
yP=[lE,fliplr(uE)];
xP=[x,fliplr(x)];
%remove nans otherwise patch won't work
xP(isnan(yP))=[];
yP(isnan(yP))=[];
H.patch=patch(xP,yP,1,'facecolor',patchColor,...
'edgecolor','none',...
'facealpha',faceAlpha);
%Make pretty edges around the patch.
H.edge(1)=plot(x,lE,'-','color',edgeColor);
H.edge(2)=plot(x,uE,'-','color',edgeColor);
%Now replace the line (this avoids having to bugger about with z coordinates)
uistack(H.mainLine,'top')
if ~holdStatus, hold off, end
if nargout==1
varargout{1}=H;
end
end
function prc = percentile(x,k)
% percentile function to replace prctile in statistics toolbox
% x .. data vector
% k .. percentage in % (k >= 1)
% if k is outside the range the min or max value of x gets assigned
len = length(x);
if isempty(x)
prc = NaN;
return
end
if len == 1
prc = x; return
end
y = sort(x);
z = 100*(0.5:1:(len-0.5))/len;
if k<z(1)
prc=k(1); return
end
if k>z(end)
prc=z(end); return
end
if isempty(find(z==k, 1))
prc = interp1(z,y,k);
else
prc = y(find(z==k, 1));
end
end
function y = nan_mean(x,dim)
% FORMAT: Y = NANMEAN(X,DIM)
%
% Average or mean value ignoring NaNs
%
% This function enhances the functionality of NANMEAN as distributed in
% the MATLAB Statistics Toolbox and is meant as a replacement (hence the
% identical name).
%
% NANMEAN(X,DIM) calculates the mean along any dimension of the N-D
% array X ignoring NaNs. If DIM is omitted NANMEAN averages along the
% first non-singleton dimension of X.
%
% Similar replacements exist for NANSTD, NANMEDIAN, NANMIN, NANMAX, and
% NANSUM which are all part of the NaN-suite.
%
% See also MEAN
if isempty(x)
y = NaN;
return
end
if nargin < 2
dim = min(find(size(x)~=1));
if isempty(dim)
dim = 1;
end
end
% Replace NaNs with zeros.
nans = isnan(x);
x(isnan(x)) = 0;
% denominator
count = size(x,dim) - sum(nans,dim);
% Protect against a all NaNs in one dimension
i = find(count==0);
count(i) = ones(size(i));
y = sum(x,dim)./count;
y(i) = i + NaN;
end
function y = nan_std(x,dim,flag)
% FORMAT: Y = NANSTD(X,DIM,FLAG)
%
% Standard deviation ignoring NaNs
%
% This function enhances the functionality of NANSTD as distributed in
% the MATLAB Statistics Toolbox and is meant as a replacement (hence the
% identical name).
%
% NANSTD(X,DIM) calculates the standard deviation along any dimension of
% the N-D array X ignoring NaNs.
%
% NANSTD(X,DIM,0) normalizes by (N-1) where N is SIZE(X,DIM). This make
% NANSTD(X,DIM).^2 the best unbiased estimate of the variance if X is
% a sample of a normal distribution. If omitted FLAG is set to zero.
%
% NANSTD(X,DIM,1) normalizes by N and produces the square root of the
% second moment of the sample about the mean.
%
% If DIM is omitted NANSTD calculates the standard deviation along first
% non-singleton dimension of X.
%
% Similar replacements exist for NANMEAN, NANMEDIAN, NANMIN, NANMAX, and
% NANSUM which are all part of the NaN-suite.
%
% See also STD
if isempty(x)
y = NaN;
return
end
if nargin < 3
flag = 0;
end
if nargin < 2
dim = min(find(size(x)~=1));
if isempty(dim)
dim = 1;
end
end
% Find NaNs in x and nanmean(x)
nans = isnan(x);
avg = nan_mean(x,dim);
% create array indicating number of element
% of x in dimension DIM (needed for subtraction of mean)
tile = ones(1,max(ndims(x),dim));
tile(dim) = size(x,dim);
% remove mean
x = x - repmat(avg,tile);
count = size(x,dim) - sum(nans,dim);
% Replace NaNs with zeros.
x(isnan(x)) = 0;
% Protect against a all NaNs in one dimension
i = find(count==0);
if flag == 0
y = sqrt(sum(x.*x,dim)./max(count-1,1));
else
y = sqrt(sum(x.*x,dim)./max(count,1));
end
y(i) = i + NaN;
end
Build with Bootstrap