The Auditory Modeling Toolbox
Go to function
LOCALIZATIONERROR - Compute psychoacoustic performance parameters for sound localization experiments
Program code:
function [varargout] = localizationerror(m,varargin)
%LOCALIZATIONERROR Compute psychoacoustic performance parameters for sound localization experiments
% Usage: [accL, precL, accP, precP, querr] = localizationerror(m)
% [res, meta, par] = localizationerror(m,errorflag)
%
% Input parameters:
% m : item list from a localization experiment. Columns:
%
% 1:4 ... azi_target,ele_target,azi_response,ele_response
%
% 5:8 ... lat_target,pol_target,lat_response,pol_response
%
% 9 ... F/B-C resolved pol_response
%
%
% LOCALIZATIONERROR(m,errorflag) returns psychoacoustic performance
% parameters for experimental response patterns.
% m is a matrix. In each row of m, the information about the target
% and response in the columns must be provided.
%
% [res, meta, par] = LOCALIZATIONERROR(m,errorflag) calculates error
% metric res given by the string errorflag. meta contains additional
% metadata, par contains additional calculation results.
%
% The errorflag may be one of:
%
% accL accuracy (average) in the lateral dimension
%
% precL precision (std. dev.) in the lateral dimension
%
% precLcentral precL derived from responses close to horizontal
% meridian (+-30 deg elevation)
%
% accP accuracy (average) in the polar dimension
%
% precP precision (circular std. dev.) in the polar dimension
%
% querr percentage of weighted polar errors > 45deg:
% [querr number_of_corrects total_number] = CalcAccPrec...
%
% accE accuracy in the elevation
%
% absaccE absolute accuracy in the elevation
%
% absaccL absolute lateral accuracy
%
% accabsL mean absolute lateral error
%
% accabsP mean absolute polar error
%
% accPnoquerr accP with quadrant errors removed
%
% precE precision in the elevation
%
% querr90 precentage of weighted polar errors > 90deg
%
% precPmedian unweighted polar precision error considering targets around
% the median plane only (+/-30deg lateral)
%
% precPmedianlocal basing on precPmedian, quadrant errors (polar error >90 deg)
% are excluded. Similar to RMS local polar error (Middlebrooks, 1999)
%
% precPnoquerr precP with quadrant errors removed
%
% rmsL lateral RMS error according to Middlebrooks (1999)
%
% rmsPmedianlocal unweighted polar RMS considering targets around
% the median plane only (+/-30deg lateral) and excluding quadrant
% errors (polar error > 90deg). Identical to RMS local polar
% error in Middlebrooks (1999).
%
% rmsPmedian unweighted polar RMS considering targets around the median plane
% only (+/-30deg lateral). Includes quadrant errors.
%
% querrMiddlebrooks quadrant errors as in Middlebrooks (1999). Use
% central positions within +/-30deg, querr is when polar
% error is >90deg
% [querr number_of_corrects total_number] = CalcAccPrec...
%
% corrcoefL lateral correlation coeffficient. Output parameter: [cc, p]
%
% corrcoefP polar correlation coeffficient, Only data within lateral
% +/-30deg are considered. Output parameter: [corr_coeff, p_of_significant_corr]
%
% SCC spherical/spatial correlation coefficient
% (Carlile et al., 1997)
%
% gainLstats performs linear regression to obtain lateral gain.
% See help of regress for a detailed description of
% the structure fields.
%
% gainL lateral gain only
%
% pVeridicalL Proportion of quasi-verdical responses in lateral
% dimension
%
% precLregress lateral scatter around linear regression line
% (only quasi-veridical responses included).
%
% sirpMacpherson2000 performs an ad-hoc selective, iterative
% regression procedure (SIRP) in order to exclude outliers and
% reversals and isolate the main concentration of responses
% in the computation of the linear fits. Outlier distance
% criterion: 40deg. Parameters [f,r] correspond to regression
% results for frontal and rear hemisphere, respectively - see
% help of regress for a detailed description of the
% structure fields.
%
% gainPfront frontal polar gain only
%
% gainPrear rear polar gain only
%
% gainP polar gain averaged between front and back
%
% slopePfront slope in degrees of regression line (frontal only)
%
% slopePrear slope in degrees of regression line (rear only)
%
% slopeP slope in degrees of regression line (front and back)
%
% pVeridicalPfront Proportion of quasi-verdical polar responses in
% the front
%
% pVeridicalPrear Proportion of quasi-verdical polar responses in
% the back
%
% pVeridicalP Proportion of quasi-verdical polar responses in
% total
%
% precPregressFront polar scatter around linear regression line for
% the front (only quasi-veridical responses included)
%
% precPregressRear same as precPregressFront but for the back
%
% precPregress average between precPregressFront and precPregressRear*
%
% perMacpherson2003 polar error rate used in Macpherson & Middlebrooks
% (2003). They measured the deviation of responses from the
% linear predictors obtained by an ad-hoc SIRP. Polar errors
% are defined by showing a deviation of >45deg
% (i.e., not being quasi-veridical) with respect
% to the linear flat stimulus prediction. Note that for this
% analysis the results from sirpMacpherson2000 are
% required and handled as LOCALIZATIONERROR(m,f,r,'perMacpherson2003')
%
% If no errorflag is provided, the function returns: accL, precL, precP, and querr
%
% Url: http://amtoolbox.sourceforge.net/data/amt-test/htdocs/amt-0.10.0/doc/general/localizationerror.php
% Copyright (C) 2009-2020 Piotr Majdak and the AMT team.
% This file is part of Auditory Modeling Toolbox (AMT) version 0.10.0
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
% References: middlebrooks1999nonindividualized macpherson2003ripples carlile1997errors
% AUTHOR : Piotr Majdak, Robert Baumgartner
%% Check input options
definput.import={'localizationerror'};
[flags,kv]=ltfatarghelper({'f','r'},definput,varargin);
if size(m,2)==1, m=m'; end
if isempty(flags.errorflag)
nargout=6;
thr=45;
accL=mean(m(:,7)-m(:,5));
varargout{1}=accL;
precL=sqrt(sum((m(:,7)-m(:,5)-accL).^2)/size(m,1));
varargout{2}=precL;
w=0.5*cos(deg2rad(m(:,7))*2)+0.5;
varargout{3}=mynpi2pi(circmean(m(:,8)-m(:,6),w));
[x,mag]=circmean(m(:,8)-m(:,6));
precP=rad2deg(acos(mag)*2)/2;
varargout{4}=precP;
idx=find((abs(mynpi2pi(m(:,8)-m(:,6))).*w)<=thr);
querr=100-size(idx,1)/size(m,1)*100;
varargout{5}=querr;
else
% if iscell(errorflag), errorflag=errorflag{1}; end
% if ~ischar(errorflag)
% error('ERR must be a string');
% end
% if size(errorflag,1)~=1, errorflag=errorflag'; end
nargout=3;
metadata=[];
par=[];
switch flags.errorflag
case 'accL'
accL=mean(m(:,7)-m(:,5));
varargout{1}=accL;
meta.ylabel='Lateral bias (deg)';
case 'absaccL'
accL=mean(m(:,7)-m(:,5));
varargout{1}=abs(accL);
meta.ylabel='Lateral absolute bias (deg)';
case 'accabsL'
accL=mean(abs(m(:,7)-m(:,5)));
varargout{1}=accL;
meta.ylabel='Lateral absolute error (deg)';
case 'corrcoefL'
[r,p] = corrcoef([m(:,5) m(:,7)]);
varargout{1}=r(1,2);
par.par1=p(1,2);
meta.ylabel='Lateral correlation coefficient';
meta.par1 = 'p for testing the hypothesis of no correlation';
case 'corrcoefP'
idx=find(abs(m(:,7))<=30);
[r,p] = corrcoef([m(idx,6) m(idx,8)]);
varargout{1}=r(1,2);
par.par1=p(1,2);
meta.ylabel='Polar correlation coefficient';
meta.par1 = 'p for testing the hypothesis of no correlation';
case 'precL'
accL=mean(m(:,7)-m(:,5));
precL=sqrt(sum((m(:,7)-m(:,5)-accL).^2)/size(m,1));
varargout{1}=precL;
meta.ylabel='Lateral precision error (deg)';
case 'precLcentral'
idx=find(abs(m(:,4))<=30); % responses close to horizontal meridian
accL=mean(m(idx,7)-m(idx,5));
precL=sqrt(sum((m(idx,7)-m(idx,5)-accL).^2)/length(idx));
varargout{1}=precL;
meta.ylabel='Central lateral precision error (deg)';
case 'rmsL'
idx=find(abs(m(:,7))<=60);
m=m(idx,:);
rmsL=sqrt(sum((mynpi2pi(m(:,7)-m(:,5))).^2)/size(m,1));
varargout{1}=rmsL;
meta.ylabel='Lateral RMS error (deg)';
case 'accP'
varargout{1}=mynpi2pi(circmean(m(:,8)-m(:,6)));
meta.ylabel='Polar bias (deg)';
case 'accabsP'
varargout{1}=mynpi2pi(circmean(abs(m(:,8)-m(:,6))));
meta.ylabel='Polar absolute bias (deg)';
case 'accPnoquerr'
w=0.5*cos(deg2rad(m(:,7))*2)+0.5;
idx=find(w>0.5);
m=m(idx,:); w=w(idx,:);
idx=find((abs(mynpi2pi(m(:,8)-m(:,6))).*w)<=45);
varargout{1}=mynpi2pi(circmean(m(idx,8)-m(idx,6)));
par.par1=length(idx);
meta.ylabel='Local polar bias (deg)';
meta.par1='Number of considered responses';
case 'precP'
w=0.5*cos(deg2rad(m(:,7))*2)+0.5;
precP=circstd(m(:,8)-m(:,6),w);
varargout{1}=precP;
meta.ylabel='Polar precision (deg)';
case 'precPnoquerr'
% w=0.5*cos(deg2rad(m(:,7))*2)+0.5;
% idx=find(w>0.5);
% thr=45; % threshold for the choice: quadrant error or not
% m=m(idx,:); w=w(idx,:);
range=360; % polar range of targets
thr=45; % threshold for the choice: quadrant error or not
latmax=0.5*180*(acos(4*thr/range-1))/pi;
idx=find(abs(m(:,7))<latmax); % targets outside not considered
m=m(idx,:);
w=0.5*cos((m(:,7))*2*pi/180)+0.5;
idx=find((abs(mynpi2pi(m(:,8)-m(:,6))).*w)<=thr);
m=m(idx,:);
w=0.5*cos((m(:,7))*pi/180*2)+0.5;
precP=circstd(m(:,8)-m(:,6),w);
varargout{1}=real(precP);
meta.ylabel='Local polar precision error (deg)';
case 'querr'
range=360; % polar range of targets
thr=45; % threshold for the choice: quadrant error or not
latmax=0.5*180*(acos(4*thr/range-1))/pi;
idx=find(abs(m(:,7))<latmax); % targets outside not considered
m=m(idx,:);
w=0.5*cos(pi*(m(:,7))/180*2)+0.5; % weigthing of the errors
corrects=size(find((abs(mynpi2pi(m(:,8)-m(:,6))).*w)<=thr),1);
totals=size(m,1);
chancerate=2*thr/range./w;
wrongrate=100-mean(chancerate)*100;
querr=100-corrects/totals*100;
varargout{1}=querr;
meta.ylabel='Weighted quadrant errors (%)';
par.par1=corrects;
meta.par1='Number of correct responses';
par.par2=totals;
meta.par2='Number of responses within the lateral range';
par.par3=wrongrate; % valid only for targets between 0 and 180°
meta.par3='Wrong rate, valid only for targets between 0 and 180°';
case 'querrMiddlebrooks' % as in Middlebrooks (199); comparable to COC from Best et al. (2005)
idx=find(abs(m(:,7))<=30);
m=m(idx,:);
idx=find((abs(mynpi2pi(m(:,8)-m(:,6))))<=90);
querr=100-size(idx,1)/size(m,1)*100;
varargout{1}=querr;
meta.ylabel='Quadrant errors (%)';
par.par1=size(idx,1);
meta.par1='Number of confusions';
par.par2=size(m,1);
meta.par2='Number of responses within the lateral range';
case 'querrMiddlebrooksRAU' % as in Middlebrooks (1999) but calculated in rau
idx=find(abs(m(:,7))<=30);
m=m(idx,:);
idx=find((abs(mynpi2pi(m(:,8)-m(:,6))))<=90);
querr=100-size(idx,1)/size(m,1)*100;
varargout{1}=rau(querr,size(m,1));
meta.ylabel='RAU quadrant errors';
par.par1=size(idx,1);
meta.par1='Number of confusions';
par.par2=size(m,1);
meta.par2='Number of responses within the lateral range';
case 'accE'
varargout{1}=mynpi2pi(circmean(m(:,4)-m(:,2)));
meta.ylabel='Elevation bias (deg)';
case 'absaccE'
varargout{1}=abs(mynpi2pi(circmean(m(:,4)-m(:,2))));
meta.ylabel='Absolute elevation bias (deg)';
case 'precE'
precE=circstd(m(:,4)-m(:,2));
varargout{1}=precE;
meta.ylabel='Elevation precision error (deg)';
case 'lateral' % obsolete
nargout=2;
varargout{1}=accL;
varargout{2}=precL;
case 'polar' % obsolete
nargout=3;
varargout{1}=accP;
varargout{2}=precP;
varargout{3}=querr;
case 'precPmedianlocal'
idx=find(abs(m(:,7))<=30);
m=m(idx,:);
idx=find(abs(mynpi2pi(m(:,8)-m(:,6)))<90);
if isempty(idx), error('Quadrant errors of 100%, local bias is NaN'); end
precPM=circstd(m(idx,8)-m(idx,6));
varargout{1}=precPM;
meta.ylabel='Local central precision error (deg)';
case 'precPmedian'
idx=find(abs(m(:,7))<=30);
m=m(idx,:);
precPM=circstd(m(:,8)-m(:,6));
varargout{1}=precPM;
meta.ylabel='Central precision error (deg)';
case 'rmsPmedianlocal'
idx=find(abs(m(:,7))<=30);
m=m(idx,:);
idx=find(abs(mynpi2pi(m(:,8)-m(:,6)))<90);
rmsPlocal=sqrt(sum((mynpi2pi(m(idx,8)-m(idx,6))).^2)/size(idx,1));
varargout{1}=rmsPlocal;
meta.ylabel='Local central RMS error (deg)';
case 'rmsPmedian'
idx=find(abs(m(:,7))<=30);
m=m(idx,:);
rmsP=sqrt(sum((mynpi2pi(m(:,8)-m(:,6))).^2)/size(m,1));
varargout{1}=rmsP;
meta.ylabel='Central RMS error (deg)';
case 'SCC'
trgt = m(:,[1 2]);
resp = m(:,[3 4]);
T = pol2dcos(trgt);
R = pol2dcos(resp);
sxy = T' * R;
sxx = T' * T;
syy = R' * R;
scc = det(sxy) / sqrt(det(sxx) * det(syy));
varargout{1}=scc;
meta.ylabel='Spatial correlation coefficient';
case 'gainLstats'
[l.b,l.bint,l.r,l.rint,l.stats]=regress(m(:,7),[ones(length(m),1) m(:,5)]);
varargout{1}=l;
meta.ylabel='Lateral gain';
case 'gainL'
l = localizationerror(m,'gainLstats');
varargout{1}=l.b(2);
meta.ylabel='Lateral gain';
case 'precLregress'
l = localizationerror(m,'gainLstats');
x = -90:90;
y = l.b(2)*m(:,5) + l.b(1); % linear regression
dev = m(:,7) - y; % deviation
dev = dev(abs(dev)<=45); % include only quasi-veridical responses
prec = rms(dev);
varargout{1} = prec;
meta.ylabel = 'Lateral scatter (deg)';
case 'pVeridicalL'
tol = 45; % tolerance of lateral angle error for quasi-veridical responses
l = localizationerror(m,'gainLstats');
x = -90:90;
y = l.b(2)*m(:,5) + l.b(1); % linear regression
dev = m(:,7) - y; % deviation wrapped to +-180 deg
prop = sum(abs(dev)<=tol)/length(dev);
varargout{1} = prop*100;
meta.ylabel = '% lateral quasi-veridical';
case 'sirpMacpherson2000'
delta = 40; % outlier tolerance in degrees
Nmin = 5; % minimum number of responses used for regression
% Front
mf = m( abs(m(:,7))<=30 & m(:,6)<90 ,:); % frontal central data only
idf = find(mf(:,8)<90); % indices correct forntal responses (init)
if length(idf)<Nmin
f.b = nan(1,2);
f.stats = nan(1,4);
else
old=[];
while not(isequal(idf,old))
[f.b,f.bint,f.r,f.rint,f.stats]=regress(mf(idf,8),[ones(length(idf),1) mf(idf,6)]);
y=f.b(2)*mf(:,6)+f.b(1);
old = idf;
dev = mod( (mf(:,8)-y) +180,360) -180; % deviation wrapped to +-180 deg
idf=find( abs(dev) < delta);
end
end
% Rear
mr = m( abs(m(:,7))<=30 & m(:,6)>=90 ,:); % rear central data only
idr = find(mr(:,8)>=90); % indices correct rear responses (init)
if length(idr)<Nmin
r.b = nan(1,2);
r.stats = nan(1,4);
else
old=[];
while not(isequal(idr,old))
[r.b,r.bint,r.r,r.rint,r.stats]=regress(mr(idr,8),[ones(length(idr),1) mr(idr,6)]);
y=r.b(2)*mr(:,6)+r.b(1);
old = idr;
dev = mod( (mr(:,8)-y) +180,360) -180; % deviation wrapped to +-180 deg
idr=find( abs(dev) < delta);
end
end
varargout{1}=f;
varargout{2}=r;
par='SIRP results';
case 'gainPfront'
f = localizationerror(m,'sirpMacpherson2000');
varargout{1}=f.b(2);
meta.ylabel='Frontal polar gain';
case 'gainPrear'
[tmp,r] = localizationerror(m,'sirpMacpherson2000');
varargout{1}=r.b(2);
meta.ylabel='Rear polar gain';
case 'gainP'
[f,r] = localizationerror(m,'sirpMacpherson2000');
varargout{1}=(f.b(2)+r.b(2))/2;
meta.ylabel='Polar gain';
case 'slopePfront'
g = localizationerror(m,'gainPfront');
varargout{1} = rad2deg(acos(1./sqrt(g.^2+1)));
meta.ylabel='Frontal polar regression slope';
case 'slopePrear'
g = localizationerror(m,'gainPrear');
varargout{1} = rad2deg(acos(1./sqrt(g.^2+1)));
meta.ylabel='Rear polar regression slope';
case 'slopeP'
g = localizationerror(m,'gainP');
varargout{1} = rad2deg(acos(1./sqrt(g.^2+1)));
meta.ylabel='Polar regression slope';
case 'pVeridicalPfront'
tol = 45; % tolerance of polar angle error for quasi-veridical responses
f = localizationerror(m,'sirpMacpherson2000');
mf = m( abs(m(:,7))<=30 & m(:,6)< 90 ,:); % frontal central data only
x = -90:270;
yf=f.b(2)*mf(:,6)+f.b(1); % linear prediction for flat, frontal targets
devf = mod( (mf(:,8)-yf) +180,360) -180; % deviation wrapped to +-180 deg
pVer = sum(abs(devf) <= tol) / length(devf);
varargout{1}=pVer*100;
meta.ylabel = '% quasi-veridical (frontal)';
case 'pVeridicalPrear'
tol = 45; % tolerance of polar angle error for quasi-veridical responses
[~,r] = localizationerror(m,'sirpMacpherson2000');
mr = m( abs(m(:,7))<=30 & m(:,6)>=90 ,:); % rear central data only
x = -90:270;
yr = r.b(2)*mr(:,6) + r.b(1); % linear prediction for flat, rear targets
devr = mod( (mr(:,8)-yr) +180,360) -180; % deviation wrapped to +-180 deg
pVer = sum(abs(devr) <= tol) / length(devr);
varargout{1}=pVer*100;
meta.ylabel = '% quasi-veridical (rear)';
case 'pVeridicalP'
tol = 45; % tolerance of polar angle error for quasi-veridical responses
[f,r] = localizationerror(m,'sirpMacpherson2000');
mf = m( abs(m(:,7))<=30 & m(:,6)< 90 ,:); % frontal central data
mr = m( abs(m(:,7))<=30 & m(:,6)>=90 ,:); % rear central data
x = -90:270;
yf = f.b(2)*mf(:,6) + f.b(1); % linear prediction for flat, frontal targets
yr = r.b(2)*mr(:,6) + r.b(1); % linear prediction for flat, rear targets
devf = mod( (mf(:,8)-yf) +180,360) -180; % deviation wrapped to +-180 deg
devr = mod( (mr(:,8)-yr) +180,360) -180; % deviation wrapped to +-180 deg
pVerf = sum(abs(devf) <= tol) / length(devf);
pVerr = sum(abs(devr) <= tol) / length(devr);
pVer = (pVerf + pVerr)/2;
varargout{1}=pVer*100;
meta.ylabel = '% quasi-veridical';
case 'precPregressFront'
tol = 45; % tolerance of polar angle error for quasi-veridical responses
f = localizationerror(m,'sirpMacpherson2000');
mf = m( abs(m(:,7))<=30 & m(:,6)< 90 ,:); % frontal central data only
x = -90:270;
yf = f.b(2)*mf(:,6) + f.b(1); % linear prediction for flat, frontal targets
devf = mod( (mf(:,8)-yf) +180,360) -180; % deviation wrapped to +-180 deg
dev = devf(abs(devf)<=tol); % include only quasi-veridical responses
prec = rms(dev);
varargout{1} = prec;
meta.ylabel = 'Frontal polar scatter (deg)';
case 'precPregressRear'
tol = 45; % tolerance of polar angle error for quasi-veridical responses
[~,r] = localizationerror(m,'sirpMacpherson2000');
mr = m( abs(m(:,7))<=30 & m(:,6)>=90 ,:); % rear central data only
x = -90:270;
yr = r.b(2)*mr(:,6) + r.b(1); % linear prediction for flat, rear targets
devr = mod( (mr(:,8)-yr) +180,360) -180; % deviation wrapped to +-180 deg
dev = devr(abs(devr)<=tol); % include only quasi-veridical responses
prec = rms(dev);
varargout{1} = prec;
meta.ylabel = 'Rear polar scatter (deg)';
case 'precPregress'
precF = localizationerror(m,'precPregressFront');
precR = localizationerror(m,'precPregressRear');
prec = (precF+precR)/2;
varargout{1} = prec;
meta.ylabel = 'Polar scatter (deg)';
case 'perMacpherson2003'
if isempty(kv.f) || isempty(kv.r)
amt_disp('Regression coefficients missing! Input is interpreted as baseline data!')
amt_disp('For this analysis the results from `sirpMacpherson2000()` for baseline data')
amt_disp('are required and must be handled as `localizationerror(m,f,r,...)`.')
[kv.f,kv.r] = localizationerror(m,'sirpMacpherson2000');
end
plotflag = false;
tol = 45; % tolerance of polar angle to be not counted as polar error
mf = m( abs(m(:,7))<=30 & m(:,6)< 90 ,:); % frontal central data only
mr = m( abs(m(:,7))<=30 & m(:,6)>=90 ,:); % rear central data only
x = -90:270;
yf=kv.f.b(2)*mf(:,6)+kv.f.b(1); % linear prediction for flat, frontal targets
yr=kv.r.b(2)*mr(:,6)+kv.r.b(1); % linear prediction for flat, rear targets
devf = mod( (mf(:,8)-yf) +180,360) -180; % deviation wrapped to +-180 deg
devr = mod( (mr(:,8)-yr) +180,360) -180;
f.pe = sum(abs(devf) > tol);
r.pe = sum(abs(devr) > tol);
per = (f.pe + r.pe) / size(m,1) * 100; % polar error rate
varargout{1} = per;
meta.ylabel = 'Polar error rate (%)';
if plotflag
plot(m(:,6),m(:,8),'ko');
axis equal
axis tight
hold on
plot(mf(:,6),yf,'LineWidth',2);
plot(mr(:,6),yr,'LineWidth',2);
plot(mf(:,6),yf+45,'r','LineWidth',2);
plot(mf(:,6),yf-45,'r','LineWidth',2);
plot(mr(:,6),yr+45,'r','LineWidth',2);
plot(mr(:,6),yr-45,'r','LineWidth',2);
title(['PER: ' num2str(per,2) '%'])
end
otherwise
error(['Unknown ERR: ' err]);
end
if length(varargout) < 2
varargout{2}=meta;
end
if length(varargout) < 3
varargout{3}=par;
end
end
end
function out_deg=mynpi2pi(ang_deg)
ang=ang_deg/180*pi;
out_rad=sign(ang).*((abs(ang)/pi)-2*ceil(((abs(ang)/pi)-1)/2))*pi;
out_deg=out_rad*180/pi;
end
function [phi,mag]=circmean(azi,dim)
% [PHI,MAG]=circmean(AZI,DIM)
%
% Calculate the average angle of AZI according
% to the circular statistics (Batschelet, 1981).
%
% For vectors, CIRCMEAN(AZI) is the mean value of the elements in AZI. For
% matrices, CIRCMEAN(AZI) is a row vector containing the mean value of
% each column. For N-D arrays, CIRCMEAN(azi) is the mean value of the
% elements along the first non-singleton dimension of AZI.
%
% circmean(AZI,DIM) takes the mean along the dimension DIM of AZI.
%
% PHI is the average angle in deg, 0..360�
% MAG is the magnitude showing the significance of PHI (kind of dispersion)
% MAG of zero shows that there is no average of AZI.
% MAG of one shows a highest possible significance
% MAG can be represented in terms of standard deviation in deg by using
% rad2deg(acos(MAG)*2)/2
%
% [PHI,MAG]=circmean(AZI,W) calculates a weighted average where each sample
% in AZI is weighted by the corresponding entry in W.
%
% Piotr Majdak, 6.3.2007
if isempty(azi) % no values
phi=NaN;
mag=NaN;
return;
end
if prod(size(azi))==1 % one value only
phi=azi;
mag=1;
return;
end
if ~exist('dim') % find the dimension
dim = min(find(size(azi)~=1));
w=ones(size(azi));
elseif prod(size(dim))~=1
% dim vector represents weights
w=dim;
dim = min(find(size(azi)~=1));
end
x=deg2rad(azi);
s=mean(sin(x).*w,dim)/mean(w);
c=mean(cos(x).*w,dim)/mean(w);
y=zeros(size(s));
for ii=1:length(s)
if s(ii)>0 & c(ii)>0
y(ii)=atan(s(ii)/c(ii));
elseif c(ii)<0
y(ii)=atan(s(ii)/c(ii))+pi;
elseif s(ii)<0 & c(ii)>0
y(ii)=atan(s(ii)/c(ii))+2*pi;
elseif s(ii)==0
y(ii)=0;
else
y(ii)=NaN;
end
end
mag=sqrt(s.*s+c.*c);
phi=wrapTo360(rad2deg(y));
end
function mag=circstd(azi,dim)
% MAG=circstd(AZI,DIM)
%
% Calculate the standard deviation angle (in deg) of AZI (in deg) according
% to the circular statistics (Batschelet, 1981).
%
% For vectors, CIRCSTD(AZI) is the std of the elements in AZI. For
% matrices, CIRCSTD(AZI) is a row vector containing the std of
% each column. For N-D arrays, CIRCSTD(azi) is the std value of the
% elements along the first non-singleton dimension of AZI.
%
% CIRCSTD(AZI,DIM) takes the std along the dimension DIM of AZI.
%
%
% MAG=CIRCSTD(AZI,W) calculates a weighted std where each vector strength
% of AZI is weighted by the corresponding entry in W.
%
% Piotr Majdak, 21.08.2008
if isempty(azi) % no values
mag=NaN;
return;
end
if prod(size(azi))==1 % one value only
mag=0;
return;
end
if ~exist('dim') % find the dimension
[phi,mag]=circmean(azi);
elseif prod(size(dim))~=1
[phi,mag]=circmean(azi,dim);
end
mag=rad2deg(sqrt(2*(1-mag)));
end
Build with Bootstrap