nt_cca_mm3

[D,E,R]=nt_cca_mm(x,y,ssize,ldaflag,nccs) - calculate metrics for match-mismatch task

  D: d-prime 
  E: error rate
  R: correlation coefficient over entire trial
  EXTRA: struct with additional info:
    .D_match: matrix of euclidean distances for matching segments/trials
    .D_mismatch: matrix of euclidean distances averaged over mismatched segments/trials
    .rms_eeg: matrix of RMS of EEG segments
    .rms_stim: matrix of RMS of stim segments

  x,y: data as trial arrays
  ssize: samples, segment size [default: all]
  ldaflag: 0: first CC, 1: LDA, 2: sum, p<1: weight with exp p
  nccs: number of CCs to keep [default: all]

0001 function [D,E,R,EXTRA]=nt_cca_mm(x,y,ssize,ldaflag,nccs)
0002 %[D,E,R]=nt_cca_mm(x,y,ssize,ldaflag,nccs) - calculate metrics for match-mismatch task
0003 %
0004 %  D: d-prime
0005 %  E: error rate
0006 %  R: correlation coefficient over entire trial
0007 %  EXTRA: struct with additional info:
0008 %    .D_match: matrix of euclidean distances for matching segments/trials
0009 %    .D_mismatch: matrix of euclidean distances averaged over mismatched segments/trials
0010 %    .rms_eeg: matrix of RMS of EEG segments
0011 %    .rms_stim: matrix of RMS of stim segments
0012 %
0013 %  x,y: data as trial arrays
0014 %  ssize: samples, segment size [default: all]
0015 %  ldaflag: 0: first CC, 1: LDA, 2: sum, p<1: weight with exp p
0016 %  nccs: number of CCs to keep [default: all]
0017 
0018 zscore=1;
0019 lda=1;
0020 
0021 if nargin<2; error('!'); end
0022 if nargin<3; ssize=[]; end
0023 if nargin<4||isempty(ldaflag); ldaflag=1; end
0024 if nargin<5; nccs=[]; end
0025 
0026 if ssize ~= round(ssize); error('!'); end
0027 
0028 % find smallest trial size across trials
0029 nsamples=size(x{1},1); % min size over trials
0030 ntrials=numel(x);
0031 for iTrial=1:ntrials
0032     if size(x{iTrial}) ~= size(y{iTrial}); error('!'); end
0033     nsamples=min(nsamples,size(x{iTrial},1));
0034 end
0035 if isempty(ssize); ssize=nsamples; end % default: segment size = nsamples
0036 nsamples=ssize*floor(nsamples/ssize); % reduce nsamples to multiple of wsize
0037 
0038 % clip all trials to smallest size
0039 if nsamples<1; error('!'); end
0040 for iTrial=1:ntrials
0041     x{iTrial}=nt_demean(x{iTrial}(1:nsamples,:)); % clip trials to new length
0042     y{iTrial}=nt_demean(y{iTrial}(1:nsamples,:));
0043 end
0044 nsegments=nsamples/ssize; % number of segments per trial
0045 
0046 % sanity check
0047 if 0 
0048     % scramble (should yield approx d-prime==0 and error == 50%)
0049     for iTrial=1:ntrials
0050         y{iTrial}=y{1+mod(iTrial+5,ntrials)};
0051         %disp([iTrial, 1+mod(iTrial+5,ntrials)]);
0052     end
0053 end
0054 
0055 % call standard nt_cca_mm to get map of good/bad segments
0056 [~,~,~,EXTRA]=nt_cca_mm(x,y,ssize,ldaflag,nccs);
0057 
0058 % determine threshold to count segment as bad
0059 d=EXTRA.DD_mismatch(:)-EXTRA.DD_match(:);  % small --> bad
0060 d=sort(d);
0061 PROP=0.1; % proportion of segments to count as bad
0062 thresh=d(round(PROP*numel(d)));
0063 
0064 DD_match=reshape(EXTRA.DD_match, [size(EXTRA.DD_match,1)/ntrials, ntrials]);
0065 DD_mismatch=reshape(EXTRA.DD_mismatch, [size(EXTRA.DD_mismatch,1)/ntrials, ntrials]);
0066 for iTrial=1:numel(x)
0067     
0068     % bad segments in this trial
0069     iBad=find(DD_mismatch(:,iTrial)-DD_match(:,iTrial)<thresh);
0070     
0071     disp(iBad);
0072     
0073     % remove bad segments
0074     xx=reshape(x{iTrial}, [ssize,size(x{iTrial},1)/ssize, size(x{iTrial},2)]); 
0075     xx(:,iBad,:)=[];
0076     xx=reshape(xx,[size(xx,1)*size(xx,2),size(xx,3)]);
0077     x2{iTrial}=xx;
0078     
0079     yy=reshape(y{iTrial},[ssize,size(y{iTrial},1)/ssize, size(y{iTrial},2)]);
0080     yy(:,iBad,:)=[];
0081     yy=reshape(yy,[size(yy,1)*size(yy,2),size(yy,3)]);
0082     y2{iTrial}=yy;
0083 end
0084 
0085 % cross-validated CCA
0086 shifts=[0]; 
0087 [AA,BB,~]=nt_cca_crossvalidate(x2,y2,shifts);
0088 
0089 % select reduced number of CCs
0090 if isempty(nccs)
0091     nccs=size(AA{1},2);
0092 else
0093     nccs=min(nccs,size(AA{1},2));
0094     for iTrial=1:ntrials
0095         AA{iTrial}=AA{iTrial}(:,1:nccs);
0096         BB{iTrial}=BB{iTrial}(:,1:nccs);
0097     end
0098     %RR=RR(1:nccs,:,:);
0099 end
0100 
0101 DD_match=[]; % normalized Euclidean distance for matching trials
0102 DD_mismatch=[]; % normalized Euclidean distance for mismatched trials
0103 RR=[];
0104 rms_eeg=[];
0105 rms_stim=[];
0106 
0107 % crossvalidation loop
0108 for iTrial=1:ntrials
0109     
0110     %{
0111     The CCA solution (AA, BB) was calculated on the basis of other trials.     
0112     We apply it to segments of this trial.
0113     %}
0114 
0115     % CCs
0116     cc_x=nt_mmat(x{iTrial},AA{iTrial});
0117     cc_y=nt_mmat(y{iTrial},BB{iTrial});
0118     %cc_x=nt_mmat(y{1+mod(iTrial,ntrials)},BB{iTrial}); % scramble
0119     
0120     % cut CCs into segments
0121     S_x=zeros(ssize,nccs,nsegments);
0122     S_y=zeros(ssize,nccs,nsegments);
0123     for iSegment=1:nsegments
0124         start=(iSegment-1)*ssize;
0125         S_x(:,:,iSegment)=cc_x(start+(1:ssize),:); 
0126         S_y(:,:,iSegment)=cc_y(start+(1:ssize),:);
0127      end
0128 
0129     % RMS of EEG and stim over each segment
0130     a=sqrt(mean(S_x.^2,1)); % average over samples
0131     a=permute(a,[2 3 4 1]); % get rid of singleton
0132     a=a(:,:); % size [nccs, nsegments*ntrials])
0133     a=mean(a); % average over CCs
0134     b=sqrt(mean(S_y.^2,1));
0135     b=permute(b,[2 3 4 1]);
0136     b=b(:,:);
0137     b=mean(b);
0138     rms_eeg=[rms_eeg;a(:)];
0139     rms_stim=[rms_stim;b(:)];
0140 
0141     
0142     % z-score each segment
0143     for iSegment=1:nsegments
0144         if zscore
0145             S_x(:,:,iSegment)=nt_normcol(nt_demean(S_x(:,:,iSegment))); % mean 0 norm 1
0146             S_y(:,:,iSegment)=nt_normcol(nt_demean(S_y(:,:,iSegment)));
0147         end
0148     end
0149     
0150     % For each envelope segment, we calculate the Euclidean distance between this segment and the corresponding
0151     % EEG segment (match) and between this segment and all the envelope segments of other trials (mismatch).
0152 
0153     % matching segments:
0154     
0155     % distance from each segment of audio to corresponding segment of EEG
0156     D_match=sqrt(mean((S_x-S_y).^2, 1));
0157     D_match=permute(D_match,[2 3 4 1]); % get rid of initial singleton
0158     D_match=D_match(:,:)'; % --> segments X comps
0159     
0160         
0161     % mismatched segments:
0162     
0163     % distance from each segment of x to all mismatched segments of y,
0164     % averaged over those segments
0165     S_x=S_x(:,:,:); 
0166     S_y=S_y(:,:,:);
0167     D_mismatch=zeros(size(S_x,3),size(S_x,3)-1, nccs);
0168     for iSegment=1:size(S_x,3)
0169         other_segments=setdiff(1:size(S_x,3),iSegment);
0170         if 1
0171             tmp=bsxfun(@minus,S_y(:,:,other_segments),S_x(:,:,iSegment));
0172         else
0173             % flip: each segment of y with all mismatched segments of x
0174             tmp=bsxfun(@minus,S_x(:,:,other_segments),S_y(:,:,iSegment));
0175         end
0176         d=sqrt(mean(tmp.^2, 1));
0177         D_mismatch(iSegment,:,:)=permute(d,[1 3 2]);
0178     end
0179     D_mismatch=mean(D_mismatch,2); % average over all other segments
0180     D_mismatch=permute(D_mismatch,[1 3 2]); 
0181     
0182     
0183     if ldaflag==1
0184         %{
0185         We want to transform distance scores (one per CC) using LDA. 
0186         
0187         To get the LDA matrix to apply to this trial, we calculate a CCA
0188         solution based on the other trials, and calculate the LDA
0189         solution from that.
0190         %}
0191          
0192         % exclude this trial
0193         other_trials=setdiff(1:ntrials,iTrial);
0194 
0195         % CCs
0196         cc_x2=nt_mmat(x(other_trials),AA{iTrial}); % AA, BB: CCA xforms derived from other trials
0197         cc_y2=nt_mmat(y(other_trials),BB{iTrial});
0198 
0199         % cut CCs into segments
0200         % !!! why not just reshape?
0201         S_x=zeros(ssize,nccs,ntrials-1,nsegments);
0202         S_y=zeros(ssize,nccs,ntrials-1,nsegments);
0203         for iTrial2=1:ntrials-1
0204             for iSegment=1:nsegments
0205                 start=(iSegment-1)*ssize;
0206 
0207                 if zscore
0208                     S_x(:,:,iTrial2,iSegment)=nt_normcol(nt_demean(cc_x2{iTrial2}(start+(1:ssize),:))); % mean 0 norm 1
0209                     S_y(:,:,iTrial2,iSegment)=nt_normcol(nt_demean(cc_y2{iTrial2}(start+(1:ssize),:)));
0210                 else
0211                     S_x(:,:,iTrial2,iSegment)=cc_x2{iTrial2}(start+(1:ssize),:); 
0212                     S_y(:,:,iTrial2,iSegment)=cc_y2{iTrial2}(start+(1:ssize),:);
0213                 end
0214             end
0215         end    
0216         % S_x and S_y are of size [nsamples, nccs, (ntrials-1), nsegments]
0217 
0218         % --> [nsamples, nccs, (ntrials-1)*nsegments]:
0219         S_x=S_x(:,:,:); 
0220         S_y=S_y(:,:,:);
0221 
0222         % For each EEG segment, we calculate the Euclidean distance between this segment and the corresponding
0223         % envelope segment (match) and between this segment and all the envelope segments of other trials (mismatch).
0224 
0225         % match:
0226         % distance from each segment of audio to corresponding segment of EEG
0227         D_match2=sqrt(mean((S_x-S_y).^2, 1));
0228         D_match2=permute(D_match2,[2 3 1]); % gets rid of initial singleton
0229         D_match2=D_match2(:,:)'; % size [nccs, (ntrials-1)*nsegments]
0230 
0231         % mismatch:
0232         % distance from each segment of audio to all mismatched segments of EEG, averaged
0233         D_mismatch2=zeros(size(S_x,3),size(S_x,3)-1, nccs);
0234         for iSegment=1:size(S_x,3)
0235             % for each segment of audio
0236             other_segments=setdiff(1:size(S_x,3),iSegment);
0237             if 0
0238                 tmp=bsxfun(@minus,S_y(:,:,other_segments),S_x(:,:,iSegment)); % sample-to-sample differences
0239             else
0240                 tmp=bsxfun(@minus,S_x(:,:,other_segments),S_y(:,:,iSegment));
0241             end
0242             d=sqrt(mean(tmp.^2, 1)); % sum squares
0243             %disp(size(tmp))
0244             D_mismatch2(iSegment,:,:)=permute(d,[1 3 2]);
0245         end
0246         D_mismatch2=mean(D_mismatch2,2); % average over all other segments
0247         D_mismatch2=permute(D_mismatch2,[1 3 2]); 
0248 
0249         if 0
0250             figure(1); clf
0251             for k=1:4
0252                 subplot (2,2,k);
0253                 histogram(D_mismatch(:,k)-D_match(:,k), -.5:.01:.5); title(mean(D_mismatch(:,k)-D_match(:,k))/std(D_mismatch(:,k)-D_match(:,k)));
0254             end
0255         end
0256 
0257         % Use DSS to find optimal linear combination of distance scores to
0258         % separate match & mismatch (--> LDA)
0259 
0260         c0=nt_cov(D_match2)/size(D_mismatch2,1);
0261         c1=nt_cov(D_mismatch2)/size(D_match2,1);
0262         [todss,pwr0,pwr1]=nt_dss0(c0,c1);
0263         if mean(D_match2*todss(:,1), 1)<0; todss=-todss; end
0264         
0265         lda_xform=todss;
0266         
0267     end % if ldaflag=1
0268 
0269     if ldaflag>0 && ldaflag<1
0270         p=ldaflag; % exponent to apply to CC distances
0271         ldaflag=3;
0272     end
0273     
0274     switch ldaflag  
0275         case 0
0276             D_match=D_match(:,1);
0277             D_mismatch=D_mismatch(:,1);
0278         case 1
0279             D_match=D_match*lda_xform(:,1);
0280             D_mismatch=D_mismatch*lda_xform(:,1);
0281         case 2
0282             D_match=mean(D_match,2);
0283             D_mismatch=mean(D_mismatch,2);
0284         case 3
0285             pp=p.^(0:size(D_match,2)-1);
0286             D_match=mean(D_match.*pp,2);
0287             D_mismatch=mean(D_mismatch.*pp,2);
0288         otherwise
0289             error('!');
0290     end
0291                 
0292     DD_match=[DD_match; D_match(:)];
0293     DD_mismatch=[DD_mismatch; D_mismatch(:)];
0294         
0295     RR(iTrial,:)=diag(corr(cc_x,cc_y));
0296 end
0297 
0298 if 0
0299     figure(100);  clf;
0300     histogram(DD_mismatch-DD_match, -.5:.05:.5); title(mean(DD_mismatch-DD_match)/std(DD_mismatch-DD_match)); 
0301     drawnow;
0302 end
0303 
0304 D=mean(DD_mismatch-DD_match, 1)/std(DD_mismatch-DD_match);
0305 E=mean(DD_mismatch-DD_match < 0, 1);
0306 R=mean(RR, 1);
0307 EXTRA.DD_mismatch=DD_mismatch;
0308 EXTRA.DD_match=DD_match;
0309 EXTRA.rms_eeg=rms_eeg;
0310 EXTRA.rms_stim=rms_stim;
0311 
0312 
0313