Home > NoiseTools > nt_peaky.m

# nt_peaky

## PURPOSE

[tocomps,ii]=nt_peaky(c,x,T,nSmooth) - find components that maximize peakiness

## SYNOPSIS

function [tocomps,ii]=nt_peaky(c,x,T,nSmooth)

## DESCRIPTION

    [tocomps,ii]=nt_peaky(c,x,T,nSmooth) - find components that maximize peakiness

      tocomps: analysis matrix
      ii: array of positions

      c: array of covariance matrices
      x: data (time*channels or time*channels*trials)
      T: (samples) time window for each covariance matrix (=dsratio)
      nSmooth: smoothing factor for covariance matrices

      This function searches over a series of time intervals for the time interval
      that maximizes the ratio between between first and second components of a DSS
      analysis that contrasts that time interval and the full data.

      Example:
        [tocomps]=nt_peaky(c): apply to array of covariance matrices
        [tocomps]=ny_peaky([],x,T): apply to signal, T is integration window
        [tocomps]=ny_peaky([],x,T,nSmooth): array is smoothed over nSmoothed samples

## CROSS-REFERENCE INFORMATION

This function calls:

- nt_dss0 [todss,pwr1,pwr2]=nt_dss0(c0,c1,keep1,keep2) - dss from covariance
- nt_greetings nt_greetings - display message the first time the toolbox is used
- nt_mmat y=nt_mmat(x,m) - matrix multiplication (with convolution)
- nt_pca0 [topcs,pwr]=nt_pca0(x,shifts,nkeep,threshold,w) - time-shift pca
- nt_xprod [y,ind]=nt_xprod(x,flag,dsratio,normrow_flag) - form all crossproducts

This function is called by:

## SOURCE CODE

    0001 function [tocomps,ii]=nt_peaky(c,x,T,nSmooth)
    0002 %[tocomps,ii]=nt_peaky(c,x,T,nSmooth) - find components that maximize peakiness
    0003 %
    0004 %  tocomps: analysis matrix
    0005 %  ii: array of positions
    0006 %
    0007 %  c: array of covariance matrices
    0008 %  x: data (time*channels or time*channels*trials)
    0009 %  T: (samples) time window for each covariance matrix (=dsratio)
    0010 %  nSmooth: smoothing factor for covariance matrices
    0011 %
    0012 %  This function searches over a series of time intervals for the time interval
    0013 %  that maximizes the ratio between between first and second components of a DSS
    0014 %  analysis that contrasts that time interval and the full data.
    0015 %
    0016 %  Example:
    0017 %    [tocomps]=nt_peaky(c): apply to array of covariance matrices
    0018 %    [tocomps]=ny_peaky([],x,T): apply to signal, T is integration window
    0019 %    [tocomps]=ny_peaky([],x,T,nSmooth): array is smoothed over nSmoothed samples
    0020 nt_greetings;
    0021
    0022 if nargin<1; error('!'); end
    0023
    0024 topcs=[];
    0025 if isempty(c);
    0026     % we need to calculate the series of covariance matrices
    0027     if nargin<3; error('!'); end
    0028     if nargin<4; nSmooth=1; end
    0029     topcs=nt_pca0(x); % improves condition of matrices
    0030     x=nt_mmat(x,topcs);
    0031     c=nt_xprod(x,'full',T);
    0032     if nSmooth>1;
    0033         c=filter(ones(nSmooth,1),1,c);
    0034         c=c(nSmooth:end,:,:);
    0035     end
    0036 end
    0037
    0038 if ndims(c)~=3 || size(c,2)~=size(c,3)
    0039     error('c has unexpected shape');
    0040 end
    0041 nComp=size(c,2);
    0042
    0043 %{
    0044 We iterate from 1 to nComp, updating the tocomps matrix at each step.
    0045 At each step we iterate over time intervals.  At each time interval we
    0046 apply DSS using interval / total power ratio as a bias.  We select the
    0047 interval that gives the greatest ratio.
    0048
    0049 We then recalculate the corresponding DSS solution, put the first component
    0050 into the tocomps matrix, and recurse on the remainder.
    0051
    0052 Calculations are done on covariance matrices rather than actual data.
    0053 %}
    0054
    0055 c0=squeeze(mean(c));
    0056 if ~isempty(topcs)
    0057     tocomps=topcs;
    0058 else
    0059     tocomps=eye(nComp);
    0060 end
    0061 iComp=1;
    0062 [m n o]= size(c); % c is ntime x nComp x nComp
    0063 while iComp<nComp
    0064
    0065     % find the index that gives the smallest D2/D1 ratio
    0066     ratio=zeros(1,size(c,1));
    0067     for iIter=1: size(c,1);
    0068         if size(c0,1)>150   % test to optimize speed
    0069             N=10;
    0070             [V,D]=eigs(squeeze(c(iIter,:,:)),c0,N);
    0071         else
    0072             [V,D]=eig(squeeze(c(iIter,:,:)),c0);
    0073         end
    0074         D=sort(diag(D),'descend');
    0075         if 0
    0076             ratio(iIter)=D(2)/D(1);
    0077         else
    0078             ratio(iIter)=mean(D(2:end))/D(1);
    0079         end
    0080     end
    0081     [~,idx]=min(ratio);
    0082     ii(iComp)=idx;
    0083
    0084     %plot(ratio); pause
    0085
    0086     %[cond(c0) cond(squeeze(c(idx,:,:)))]
    0087
    0088     % DSS
    0089     [todss,pwr0,pwr1]=nt_dss0(c0,squeeze(c(idx,:,:)),[],10.^-10);
    0090     profile=pwr1./pwr0;
    0091     score(iComp)=profile(1)/profile(2);
    0092
    0093     % update tocomps
    0094     tocomps=[tocomps(:,1:iComp-1), tocomps(:,iComp:end)*todss];
    0095     nComp=size(tocomps,2);
    0096
    0097     % update covariance matrices
    0098     if iComp<nComp
    0099         todss=todss(:,2:end);
    0100         c0=todss'*c0*todss;
    0101         %{
    0102         For speed, we reshape c to apply left and right multiplication to all
    0103         matrices at once.
    0104         %}
    0105         [m n o]=size(c);
    0106         c=reshape(permute(c,[2 3 1]), [n o*m]); % concatenate horizontally
    0107         c=todss'*c; % left multiply
    0108         n=size(c,1);
    0109         c=reshape(permute(reshape(c,[n o m]),[1 3 2]),[m*n o]); % concatenate vertically
    0110         c=c*todss; % right multiply
    0111         o=size(c,2);
    0112         c=permute(reshape(c, [n m o]),[2 1 3]) ; % reshape to ntime x nComp x nComp
    0113     end
    0114     iComp=iComp+1;
    0115 end
    0116
    0117
    0118
    0119