analysisPlot

 function [] = analysisPlot(rr,K,fr0,ind,name,units,rms,cla,phi)

0001 % function [] = analysisPlot(rr,K,fr0,ind,name,units,rms,cla,phi)
0002 function [] = analysisPlot(rr,K,fr0,ind,opt)
0003 
0004 par = {'name','units','rms','cla','phi','mds','fps','trng','fg'};
0005 def = {'','','[]','[]','[]','[1,3]','500','[-50,150]',num2str(ceil(20*rand))};
0006 cla = [];
0007 
0008 for k = 1:length(par)
0009   p = par{k};
0010   d = def{k};
0011   eval([p,' = [];']);
0012   if isfield(opt,p)
0013     eval([p,' = opt.',p,';']);
0014   else
0015     eval([p,' = ',d,';']);
0016   end
0017 end
0018 
0019 % keyboard
0020 
0021 if isempty(K)
0022   K = 1:length(rr);
0023 end
0024 
0025 % keyboard
0026 
0027 rr = rr(K);
0028 
0029 if nargin < 9
0030   sc = 1000/fps;
0031   vsphi = 0;
0032   phi = 0;
0033   s = ((1:size(rr(1).atd,2)) - fr0 + 1).*sc;
0034   xlim = trng.*sc;
0035   xlab = 'Time (ms)';
0036 else
0037   vsphi = 1;
0038   s = phi;
0039   xlim = pi*[-5,8];
0040   xlab = 'Phase (rad)';
0041 end
0042 s0 = find(s <= xlim(1),1,'last');
0043 sf = find(s >= xlim(2),1,'first');
0044 
0045 if isempty(cla)
0046   docla = 0;
0047   cla = 0;
0048 else
0049   docla = 1;
0050   cla = cla(K);
0051 end
0052 
0053 if nargin < 7
0054   rms = 0;
0055 if nargin < 6
0056   units = '';
0057 end
0058 end
0059 
0060 if rms
0061   name = ['RMS ',name];
0062 end
0063 
0064 if docla
0065   name = ['Cla ',name];
0066 end
0067 
0068 if vsphi
0069   name = [name,' vs Phase'];
0070 end
0071 
0072 bname = name;
0073 
0074 % s = ((1:size(rr(1).atd,2)) - fr0 + 1).*2;
0075 % xlim = [-50,150].*2;
0076 
0077 ps = {{'b','b','r','g'},{'r','c','m','y'}};
0078 % pps = {{[0,0,0.5],[0.5,0,0],[0,0.5,0]},{[1,0.5,0.5],[.5,1,.5],[.5,.5,1]}};
0079 
0080 
0081 loc = 'SouthWest';
0082 
0083 ylim = [nan,nan];
0084 
0085 % close all;
0086 for md = mds
0087   fig = figure(fg);
0088   clf;
0089   if md == 3
0090 %     set(fig,'Position',[100,250,800,300]);
0091 %     name = ['Bootstrap Residual ',bname];
0092     name = ['Residual ',bname];
0093   elseif md == 2
0094 %     set(fig,'Position',[100,50,800,300]);
0095     name = ['Residual ',bname];
0096   else
0097 %     set(fig,'Position',[100,500,800,300]);
0098   end
0099   hold on;
0100 
0101   leg = {'Original','Control','Mass','Inertia'};
0102   L = {};
0103   for t = unique([rr.tmt])
0104     for k = 1:(2^docla)
0105       L{end+1} = leg{t+1};
0106       plot(xlim(1)+[-2;-1],[0;0],ps{k}{t+1},'LineWidth',3);
0107     end
0108   end
0109   leg = legend(L,'Location',loc);
0110   
0111   blk = plot(xlim(1)+[-2;-1],[0;0],'--','LineWidth',2,'color',0.2*[1,1,1]);
0112 
0113   N = [];
0114   Na = [];
0115 
0116   for t = unique([rr.tmt])
0117     for k = 1:(2^docla)
0118       K = find([rr.tmt] == t & (~docla | (cla == k)));
0119       N(end+1) = length(K);
0120       Na(end+1) = length(unique([rr(K).id]));
0121 
0122       if isempty(K)
0123         continue;
0124       end
0125 
0126       if ~vsphi
0127         td   = cat(3,rr(K).atd);
0128         if docla
0129           mdXV = cat(3,rr(K).amdXVc);
0130         else
0131           mdXV = cat(3,rr(K).amdXV);
0132         end
0133       else
0134         td   = cat(3,rr(K).atdp);
0135         mdXV = cat(3,rr(K).amdXVp);
0136       end
0137       
0138       if any(md == [2,3])
0139         rd = td - mdXV;
0140       else
0141         rd = td;
0142       end
0143 
0144       if rms
0145         d = sqrt(squeeze(nanmean(rd(ind,:,:).^2,1)));
0146       else
0147         d = squeeze(nanmean(rd(ind,:,:),1));
0148       end
0149       
0150 %       figure;
0151 %       plot(d);
0152 %       legend(num2str((1:size(d,2))'))
0153 %       keyboard
0154       
0155       D = d;
0156       J = find(isnan(d(s0,:)) | isnan(d(sf,:)));
0157       
0158 %       disp(['xlim killed ',num2str(length(J)),' trials']);
0159       
0160 %       keyboard
0161       
0162       N(end) = length(setdiff(K,J));
0163       Na(end) = length(unique([rr(setdiff(K,J)).id]));
0164       
0165       d(:,J) = [];
0166       
0167 %       if nanmean(nanmean(d(1:fr0-20,:),2)) < nanmean(nanmean(d(fr0+20:fr0+50,:),2))
0168 %         loc = 'NorthWest';
0169 %       end
0170       
0171 %       d(s0:sf,:) = filtfilt(B,A,d(s0:sf,:));
0172       
0173 %       keyboard
0174 
0175       if 1 %md == 3
0176         
0177         dd = nan*d(:,1);
0178         for b = 1:200
0179           bi = floor(size(d,2)*rand(2*size(d,2),1))+1;
0180 
0181           dd(:,b) = nanmean(d(:,bi),2);
0182         end
0183 
0184         da = nanmean(dd,2);
0185 %         de = [prctile(dd,2.5,2),prctile(dd,97.5,2)];
0186         de = [prctile(dd,25,2),prctile(dd,75,2)];
0187         
0188         ylim = [min([ylim(1);de(s0:sf,1)]),max([ylim(2);de(s0:sf,2)])];
0189         
0190 %         keyboard
0191 
0192         plot(s, smth(da), ps{k}{t+1}, 'LineWidth', 2);
0193         plot(s, smth(de), [ps{k}{t+1},'--'], 'LineWidth',1);
0194       else
0195       
0196         av = nanmean(d,2);
0197 % %       Standard error tests whether the populations are different
0198 % %       sta = 2*nanstd(d,[],2)./sqrt(Na(end));
0199 %       st = 2*nanstd(d,[],2)./sqrt(N(end));
0200 % %       Standard deviation tests whether the means are different
0201 % %       sta = 2*nanstd(d,[],2)./Na(end);
0202 %       st = 2*nanstd(d,[],2)./N(end);
0203         pr = [prctile(d,25,2),prctile(d,75,2)];
0204 
0205 %       plot(s, d,  'color', pps{k}{t});
0206         plot(s, av, ps{k}{t+1}, 'LineWidth',2);
0207 %       plot([s;s]', [av,av]+[-sta,+sta], [ps{k}{t},':'], 'LineWidth',1);
0208 %       plot([s;s]', [av,av]+[-st,+st], [ps{k}{t},'--'], 'LineWidth',1);
0209         plot(s, pr, [ps{k}{t+1},'--'], 'LineWidth',1);
0210         
0211         ylim = [min([ylim(1);pr(s0:sf,1)]),max([ylim(2);pr(s0:sf,2)])];        
0212 
0213       end
0214     end
0215   end
0216       
0217   if nanmean(nanmean(d(1:fr0-20,:),2)) < mean(ylim)
0218     loc = 'NorthWest';
0219   end
0220 
0221   str = get(leg,'String');
0222   for k = 1:length(str)
0223     str{k} = [str{k},', N = ',num2str(N(k))];
0224   end
0225   set(leg,'String',str,'Location',loc);
0226   
0227   ylim = ylim+diff(ylim)*0.05*[-1,1];
0228 
0229 %   ylim = get(gca,'YLim');
0230   set(blk,'XData',[xlim,nan,0,0],'YData',[0,0,nan,ylim]);
0231   
0232   set(gca,'XLim',xlim,'YLim',ylim,'FontName','Times','FontSize',12);
0233   xlabel(xlab,'FontName','Times','FontSize',14);
0234   ylabel([name,' ',units],'FontName','Times','FontSize',14);
0235 
0236   grid on;
0237 
0238   set(fig,'PaperPositionMode','auto');
0239 
0240   sname = name;
0241   sname(sname == ' ') = '_';
0242 
0243   saveas(fig,['fig/',sname,'.png'],'png');
0244   saveas(fig,['fig/',sname,'.pdf'],'pdf');
0245 end
0246 
0247 
0248 function s = smth(n)
0249 
0250 [B,A] = butter(2,0.2);
0251 nn = ~isnan(n);
0252 s = n;
0253 s(nn) = filtfilt(B,A,n(nn));
0254