Home > NoiseTools > nt_pca_kmeans.m

# nt_pca_kmeans

## PURPOSE

[topcs,pwr]=nt_pca_kmeans(x,nkeep) - PCA preceded by kmeans for speed

## SYNOPSIS

function [topcs,pwr]=nt_pca_kmeans(x,shifts,nkeep)

## DESCRIPTION

    [topcs,pwr]=nt_pca_kmeans(x,nkeep) - PCA preceded by kmeans for speed

      topcs: PCA matrix

      x: data (time*channels or time*channels*trials)
      nkeep: desired number of PCs

     The kmeans implementation of VLFeat appears to be very efficient,
     so we use it to find a small set of clusters to which we apply PCA.

     This differs from the usual way of combining PCA and kmeans, which is
     to use PCA to initialize kmeans.

## CROSS-REFERENCE INFORMATION

This function calls:

- nt_mmat y=nt_mmat(x,m) - matrix multiplication (with convolution)
- nt_multishift z=nt_multishift(x,shifts,amplitudes) - apply multiple shifts to matrix
- nt_pca0 [topcs,pwr,y]=nt_pca0(x,shifts,nkeep,threshold,w) - time-shift pca
- nt_unfold y=nt_fold(x) - unfold 3D to 2D

This function is called by:

## SOURCE CODE

    0001 function [topcs,pwr]=nt_pca_kmeans(x,shifts,nkeep)
    0002 %[topcs,pwr]=nt_pca_kmeans(x,nkeep) - PCA preceded by kmeans for speed
    0003 %
    0004 %  topcs: PCA matrix
    0005 %
    0006 %  x: data (time*channels or time*channels*trials)
    0007 %  nkeep: desired number of PCs
    0008 %
    0009 % The kmeans implementation of VLFeat appears to be very efficient,
    0010 % so we use it to find a small set of clusters to which we apply PCA.
    0011 %
    0012 % This differs from the usual way of combining PCA and kmeans, which is
    0013 % to use PCA to initialize kmeans.
    0014
    0015 if nargin<3; error('!'); end
    0016 if isempty(shifts); shifts=0; end
    0017 if shifts
    0018     x=nt_multishift(x,shifts);
    0019 end
    0020 if nkeep>size(x,2); error('!'); end
    0021
    0022 x=nt_unfold(x);
    0023 [nsamples,nchans]=size(x);
    0024
    0025 if nkeep>nsamples; error('nkeep greater than number of samples'); end
    0026 if nkeep>nchans; error('nkeep greater than number of channels'); end
    0027
    0028
    0029 % normalize channels
    0030 nrm=sqrt(mean(x.^2));
    0031 tonrm=diag(1./nrm); % normalization matrix
    0032 tonrm(find(isnan(x)))=0;
    0033 x=x*tonrm;
    0034
    0035 % cluster
    0036 nclusters=round(sqrt(nchans));
    0037 [C,A]=vl_kmeans(x,nclusters,'algorithm','elkan');
    0038
    0039 % perform PCA on each cluster
    0040 topcs1=zeros(nchans);
    0041 pwr=zeros(1,nchans);
    0042 idx=0;
    0043 for k=1:nclusters
    0044     %[nclusters k numel(find(A==k))]
    0045     [m,p]=nt_pca0(x(:,find(A==k)));
    0046     topcs1(find(A==k),idx+(1:size(m,2))) = m;
    0047     pwr(idx+(1:size(m,2)))=p;
    0048     idx=idx+size(m,2);
    0049 end
    0050
    0051 % sort by decreasing power
    0052 [~,idx]=sort(pwr,'descend');
    0053
    0054 topcs2=nt_pca0(nt_mmat(x,topcs1(:,idx(1:nkeep))));
    0055
    0056 topcs=tonrm*(topcs1(:,idx(1:nkeep))*topcs2);
    0057
    0058
    0059
    0060
    0061
    0062
    0063
    0064