Home > NoiseTools > nt_kurtosis.m

# nt_kurtosis

## PURPOSE

[todss,K]=nt_kurtosis(x,nIterations,exponent,w,smooth)- find high-kurtosis components

## SYNOPSIS

function [todss,K]=nt_kurtosis(x,nIterations,exponent,w,smooth)

## DESCRIPTION

    [todss,K]=nt_kurtosis(x,nIterations,exponent,w,smooth)- find high-kurtosis components

     x: data
     nIterations: number of times to iterate calculation (default: 5)
     exponent: exponent to apply to mask [default: 1]
     w: weighting function
     smooth: samples, smoothing to apply (twice) to mask

## CROSS-REFERENCE INFORMATION

This function calls:

- nt_cov [c,tw]=nt_cov(x,shifts,w) - time shift covariance
- nt_demean [y,mn]=nt_demean(x,w) - remove weighted mean over cols
- nt_dss0 [todss,pwr1,pwr2]=nt_dss0(c0,c1,keep1,keep2) - dss from covariance
- nt_kurtosis [todss,K]=nt_kurtosis(x,nIterations,exponent,w,smooth)- find high-kurtosis components
- nt_mmat y=nt_mmat(x,m) - matrix multiplication (with convolution)
- nt_unfold y=nt_fold(x) - unfold 3D to 2D

This function is called by:

- nt_kurtosis [todss,K]=nt_kurtosis(x,nIterations,exponent,w,smooth)- find high-kurtosis components

## SOURCE CODE

    0001 function [todss,K]=nt_kurtosis(x,nIterations,exponent,w,smooth)
    0002 %[todss,K]=nt_kurtosis(x,nIterations,exponent,w,smooth)- find high-kurtosis components
    0003 %
    0004 % x: data
    0005 % nIterations: number of times to iterate calculation (default: 5)
    0006 % exponent: exponent to apply to mask [default: 1]
    0007 % w: weighting function
    0008 % smooth: samples, smoothing to apply (twice) to mask
    0009
    0010 if nargin<1; error('!'); end
    0011 if nargin<2||isempty(nIterations); nIterations=5; end
    0012 if nargin<3||isempty(exponent); exponent=2; end
    0013 if nargin<4; w=[]; end
    0014 if nargin<5; smooth=[];end
    0015
    0016 if ndims(x)>2; x=nt_unfold(x); end
    0017
    0018 KK={};
    0019 K0=kurtosis(x);
    0020 if nIterations==1
    0021
    0022     %{
    0023     We choose the channel with highest kurtosis to define a mask to
    0024     emphasize high amplitude intervals. DSS then finds components with
    0025     variance maximal within those intervals.
    0026     %}
    0027
    0028     K=kurtosis(x);
    0029     [~,idx]=max(K);
    0030     mask=x(:,idx).^exponent;
    0031     mask=filtfilt(ones(5,1),1,mask);
    0032     if ~isempty(smooth)&&smooth>1;
    0033         mask=filtfilt(ones(1,smooth),1,mask); % filtfilt keeps mask aligned with data
    0034     end
    0035     c0=nt_cov(x,[],w);
    0036     c1=nt_cov(nt_demean(x.*bsxfun(@times,x,mask),w),[],w);
    0037     [todss]=nt_dss0(c0,c1);
    0038     KK{1}=kurtosis(nt_mmat(x,todss));
    0039
    0040 else
    0041
    0042     %{
    0043     Apply DSS repeatedly, concatenating the DSS matrices. We hope (no
    0044     guarantee) that the solution will reveal components with high kurtosis.
    0045     %}
    0046
    0047     z=x;
    0048     todss=eye(size(x,2));
    0049     for k=1:nIterations
    0050         [T]=nt_kurtosis(z,1,exponent,w,smooth);
    0051         z=nt_mmat(z,T);
    0052         todss=todss*T;
    0053         KK{k}=kurtosis(z);
    0054     end
    0055
    0056 end
    0057 K=KK{end};
    0058
    0059 if nargout==0;
    0060     % don't return anything, just plot
    0061     figure(101); clf
    0062     subplot 121;
    0063     plot(K0); hold on
    0064     for k=1:nIterations
    0065         plot(KK{k});
    0066     end
    0067     nt_colorlines;
    0068     title('kurtosis score for each iteration'); xlabel('component');
    0069     legend(num2str((0:nIterations)')); set(gca,'yscale','log')
    0070     subplot 222;
    0071     z=nt_mmat(x,todss(:,1));
    0072     plot(z); title('first DSS component');
    0073     subplot 224;
    0074     z=nt_mmat(x,todss(:,end));
    0075     plot(z); title('last DSS component');
    0076     clear c0 c1 todss pwr0 pwr1 K
    0077 end
    0078