Home > NoiseTools > nt_enforce_stationarity.m

# nt_enforce_stationarity

## PURPOSE

y=nt_enforce_stationarity(x,DSR,thresh) - locally project out non-stationary components

## SYNOPSIS

function x=nt_enforce_stationarity(x,DSR,thresh);

## DESCRIPTION

     y=nt_enforce_stationarity(x,DSR,thresh) - locally project out non-stationary components

     y: cleaned data

     x: data to clean (time * channels)
     DSR: size of chunks over which to cluster
     thresh: variance ratio threshold

## CROSS-REFERENCE INFORMATION

This function calls:

- nt_cluster_jd [IDX,todss,SCORE]=nt_cluster_jd(x,dsr,smooth,flags,init,verbose) - cluster with joint diagonalization
- nt_cov [c,tw]=nt_cov(x,shifts,w) - time shift covariance
- nt_dss0 [todss,pwr1,pwr2]=nt_dss0(c0,c1,keep1,keep2) - dss from covariance
- nt_fold y=fold(x,epochsize) - fold 2D to 3D
- nt_mmat y=nt_mmat(x,m) - matrix multiplication (with convolution)
- nt_pca [z,idx]=nt_pca(x,shifts,nkeep,threshold,w) - time-shift pca
- nt_tsr_nodemean [y,idx,w]=nt_tsr_nodemean(x,ref,shifts,wx,wref,keep,thresh) - time-shift regression (TSPCA)
- nt_unfold y=nt_fold(x) - unfold 3D to 2D

This function is called by:

## SOURCE CODE

    0001 function x=nt_enforce_stationarity(x,DSR,thresh);
    0002 % y=nt_enforce_stationarity(x,DSR,thresh) - locally project out non-stationary components
    0003 %
    0004 % y: cleaned data
    0005 %
    0006 % x: data to clean (time * channels)
    0007 % DSR: size of chunks over which to cluster
    0008 % thresh: variance ratio threshold
    0009
    0010 if nargin<3 || isempty(thresh); thresh=10; end
    0011 if nargin<2 || isempty(DSR); error('!'); end
    0012
    0013 nx=size(x,1);
    0014 x=nt_unfold(x);
    0015
    0016
    0017 score=inf;
    0018
    0019 figure(10); clf; set(gcf, 'name', 'nt_enforce_stationarity');
    0020 subplot 413; plot(x); title('original')
    0021
    0022 while 1
    0023
    0024     % find clusters with maximally different covariance
    0025     [A]=nt_cluster_jd((nt_pca(x)),DSR);
    0026     figure(10);
    0027     subplot 411; plot(A,'.-'); title('cluster mask'); drawnow
    0028
    0029     % DSS to enhance the smaller cluster
    0030     if numel(find(A==1))<numel(find(A==2));
    0031         idx=find(A==1);
    0032         [todss]=nt_dss0(nt_cov(x),nt_cov(x(find(A==1),:)));
    0033     else
    0034         idx=find(A==2);
    0035         %todss=fliplr(todss);
    0036         [todss]=nt_dss0(nt_cov(x),nt_cov(x(find(A==2),:)));
    0037     end
    0038
    0039     z=nt_mmat(x,todss);
    0040     subplot 412; plot(z(:,1)); title('component to remove');
    0041     score=mean(z(idx,1).^2)/mean(z(:,1).^2);
    0042     disp(['power ratio score: ', num2str(score)]);
    0043
    0044     if score<thresh
    0045         break
    0046     end
    0047
    0048     x(idx,:)=nt_tsr_nodemean(x(idx,:),z(idx,1));
    0049     subplot 414; plot(x); title('clean');
    0050     drawnow
    0051
    0052     %pause;
    0053     %x=nt_demean(x);
    0054 end
    0055 x=nt_fold(x,nx);
    0056 if nargout==0;
    0057     clear x
    0058 end
    0059
    0060
    0061