test_nt_kurtosis

 Test nt_bias_kurtosis on synthetic data

0001 
0002 % Test nt_bias_kurtosis on synthetic data
0003 clf;
0004 rng(0); % reset random number generator
0005 nsamples=5000;
0006 nchans=10;
0007 
0008 %{
0009 We create a peaky target and mix it with multiple sources of noise via a
0010 random mixing matrix.
0011 %}
0012 
0013 target=zeros(nsamples,1);
0014 target(900:1100,:)=1; %target(1001:1100,:)=-1;
0015 noise=randn(nsamples,nchans);
0016 x=noise(:,1:nchans-1)*randn(nchans-1,nchans);
0017 SNR=0.001;
0018 x=x + SNR*target*randn(1,nchans);
0019 x=nt_demean(x);
0020 
0021 figure(1); clf
0022 subplot 221; plot (target); title('target'); ylim([0 1.1]);
0023 subplot 222; plot(x); title('mixture');
0024 
0025 %{
0026 We use nt_bias_kurtosis to find a transform that reveals components with
0027 high kurtosis.
0028 %}
0029 nIterations=5;
0030 exponent=2;
0031 w=[];
0032 smooth=1;
0033 [todss]=nt_kurtosis(x,nIterations,exponent,w,smooth);
0034 
0035 z=nt_mmat(x,todss);
0036 subplot 224
0037 plot(z(:,1)); title('extracted');
0038 
0039 %{
0040 Same, with two targets.
0041 %}
0042 
0043 target1=zeros(nsamples,1);
0044 target1(900:1000,:)=1; %target(1001:1100,:)=-1;
0045 target2=zeros(nsamples,1);
0046 target2(2001:2200,:)=1; %target(1001:1100,:)=-1;
0047 x=noise(:,1:nchans-2)*randn(nchans-2,nchans);
0048 SNR=0.01;
0049 x=x + SNR*target1*randn(1,nchans) + SNR*target2*randn(1,nchans);
0050 x=nt_demean(x);
0051 
0052 figure(2); clf
0053 subplot 221; plot ([target1,target2]); title('targets'); ylim([0 1.1]);
0054 subplot 222; plot(x); title('mixture');
0055 
0056 nIterations=5;
0057 exponent=2;
0058 w=[];
0059 smooth=1;
0060 [todss]=nt_kurtosis(x,nIterations,exponent,w,smooth);
0061 z=nt_mmat(x,todss);
0062 
0063 subplot 426
0064 plot(z(:,1)); title('first peaky component');
0065 
0066 %{
0067 We use the peakiest component to define a weight that deemphasizes the
0068 interval for which that component has high amplitude.  
0069 We then call nt_bias_kurtosis to find remaining peaky components.
0070 %}
0071 
0072 w=abs(z(:,1))<0.5*max(abs(z(:,1)));
0073 nIterations=5;
0074 exponent=2;
0075 smooth=1;
0076 [todss]=nt_kurtosis(x,nIterations,exponent,w,smooth);
0077 z=nt_mmat(x,todss);
0078 
0079 subplot 428;
0080 plot(z(:,1)); title('second peaky component');
0081