test_nt_star

0001 clear
0002 
0003 %{
0004 Test 1:
0005 N channels,  1 sinusoidal target, N-3 noise sources, temporally local
0006 artifacts on each channel.
0007 %}
0008 
0009 nsamples=1000;
0010 f=2;
0011 target=sin((1:nsamples)'/nsamples*2*pi*f);
0012 nchans=10;
0013 noise=randn(nsamples,nchans-3);
0014 
0015 SNR=sqrt(1);
0016 x0=nt_normcol(noise*randn(size(noise,2),nchans))+ SNR * target*randn(1,nchans);
0017 x0=nt_demean(x0);
0018 artifact=zeros(size(x0));
0019 for k=1:nchans
0020     artifact( (k-1)*100 + (1:20) , k)=1;
0021 end
0022 x=x0+20*artifact;
0023 
0024 [y,w]=nt_star(x,2); 
0025 figure(1); clf; subplot 311; plot(x); subplot 312;  plot(y); title(['SNR= ' ,num2str(SNR)]);
0026 subplot 313; plot(nt_demean(y-x0));
0027 
0028 SNR=sqrt(10^-7);
0029 x0=nt_normcol(noise*randn(size(noise,2),nchans))+ SNR * target*randn(1,nchans);
0030 x0=nt_demean(x0);
0031 artifact=zeros(size(x0));
0032 for k=1:nchans
0033     artifact( (k-1)*100 + (1:20) , k)=1;
0034 end
0035 x=x0+20*artifact;
0036 
0037 %x=[x,x(:,1)]; % make data rank deficient
0038 
0039 [y,w]=nt_star(x,2); 
0040 
0041 figure(2); clf; subplot 311; plot(x);title(['SNR= ' ,num2str(SNR)]);
0042 subplot 312;  plot(y); title('artifact removed');
0043 subplot 313; plot(nt_demean(y-x0)); title('error after artifact removal')
0044 
0045 
0046 x=nt_demean(x);
0047 y=nt_demean(y);
0048 
0049 c0=nt_cov(x);
0050 c1=nt_cov(nt_detrend(x,2,[],'sinusoids'));
0051 todss=nt_dss0(c0,c1);
0052 z1=nt_normcol(nt_mmat(x,todss));
0053 c0=nt_cov(y);
0054 c1=nt_cov(nt_detrend(y,2,[],'sinusoids'));
0055 todss=nt_dss0(c0,c1);
0056 z2=nt_normcol(nt_mmat(y,todss));
0057 
0058 figure(3); clf; set(gcf, 'color',[1 1 1], 'position', [300   500   520 600])
0059 
0060 subplot 511
0061 plot(target)
0062 set(gca,'fontsize', 14, 'xtick',[], 'ytick', []);
0063 ylim([-1.1 1.1]); title('target')
0064 
0065 subplot 512
0066 plot(x)
0067 set(gca,'fontsize', 14, 'xtick',[], 'ytick', []);
0068 mn=min(x(:)); mx=max(x(:));
0069 ylim([mn-(mx-mn)*0.1, mx+(mx-mn)*0.1])
0070 title('mixed with noise and artifacts')
0071 
0072 subplot 513
0073 plot(y)
0074 set(gca,'fontsize', 14, 'xtick',[], 'ytick', []);
0075 title('STAR')
0076 
0077 subplot 515
0078 plot(z2(:,end));
0079 mn=min(z2(:,end)); mx=max(z2(:,end));
0080 set(gca,'fontsize', 14, 'ytick', []); xlabel('samples');
0081 ylim([mn-(mx-mn)*0.1, mx+(mx-mn)*0.1])
0082 %text(800,mx-(mx-mn)*0.1, 'SNR=10^-8','fontsize',14)
0083 title('STAR + JD')
0084 
0085 subplot 514
0086 plot(z1(:,end));
0087 set(gca,'fontsize', 14, 'xtick',[], 'ytick', []);
0088 mn=min(z1(:,end)); mx=max(z1(:,end));
0089 ylim([mn-(mx-mn)*0.1, mx+(mx-mn)*0.1])
0090 title('JD')
0091 
0092 
0093 
0094