test_nt_qca1

{

0001 %{
0002 
0003 Extract component based on reproducible induced activity
0004 
0005 %}
0006 
0007 reset(RandStream.getGlobalStream); % to get reproducible signals
0008 
0009 sr=1000;
0010 nsamples=1000;
0011 ntrials=100;
0012 nchans=10;
0013 CF1=4;
0014 CF2=2; 
0015 FSIZE=12;
0016 
0017 % sinusoidal pulse
0018 target=sin(2*pi*0.5*(1:nsamples)'/sr).^2;
0019 target=repmat(target,[1,1,ntrials]) .* sin( 2*pi*(CF1*repmat((1:nsamples)',[1,1,ntrials])/sr + repmat(rand(1,1,ntrials),[nsamples,1,1])));
0020 target=nt_normcol(target);
0021 
0022 % noise
0023 NNOISE=9;
0024 noise=randn(nsamples,NNOISE,ntrials);
0025 noise=nt_mmat(noise,randn(NNOISE,nchans));
0026 noise=nt_normcol(noise);
0027 
0028 % data
0029 SNR=0.0001;
0030 x=noise+ SNR * nt_mmat(target,randn(1,nchans)) ;
0031 
0032 x=nt_demean(x);
0033 x=nt_normcol(x);
0034 
0035 [squares,quads,D]=nt_qca(x,[],[],10);
0036 
0037 
0038 figure(1); clf;
0039 set(gca,'fontsize',12)
0040 set(gcf,'color',[1 1 1])
0041 set(gcf, 'position', [667   368   700   400])
0042 
0043 subplot 231
0044 plot(0.9*squeeze(target(:,1,1:5)/max(target(:))), 'k'); set(gca,'ytick',[]); title('target (5 trials)', 'fontsize',14);
0045 subplot 232
0046 plot(squeeze(x(:,1,:))); set(gca,'ytick',[]); title('mixture', 'fontsize',14); xlabel('samples', 'fontsize', 14);
0047 subplot 233
0048 plot(0.9*squeeze(squares(:,1,1:5)/max(squares(:))), 'k');  set(gca,'ytick',[]); title('recovered', 'fontsize',14);
0049 
0050 clear x; clear noise
0051 
0052 
0053 %return
0054 P=mfilename('fullpath');
0055 [PATHSTR,NAME,EXT] = fileparts(P);
0056 
0057 % noise is meg data
0058 load([PATHSTR,'/../DATA/meg']);
0059 meg=nt_unfold(meg);
0060 meg=nt_fold(meg(1:76000,:),1000);
0061 meg=nt_demean2(meg);
0062 sr=1000; % Hz
0063 
0064 [idx,d]=nt_find_outlier_trials2(meg,1.5);% plot(d);
0065 meg=meg(:,:,idx);
0066 [idx,d]=nt_find_outlier_trials2(meg,1.5); % plot(d);
0067 meg=meg(:,:,idx);
0068 nsamples=size(meg,1); nchans=size(meg,2); ntrials=size(meg,3);
0069 
0070 
0071 % add 'target'
0072 if 0
0073     target=sin(2*pi*0.5*(1:nsamples)'/sr).^2;
0074     target=repmat(target,[1,1,ntrials]).*sin(2*pi*(20*repmat((1:nsamples)',[1,1,ntrials])/sr+repmat(rand(1,1,ntrials),[nsamples,1,1])));
0075 else
0076     % target is pulses with random latency
0077     PWIDTH=60; PRANGE=400;
0078     B=sin(2*pi*(1:PWIDTH)/(2*PWIDTH));
0079     target=zeros(nsamples,1,ntrials);
0080     mixmatrix=randn(1,nchans);
0081     for k=1:ntrials
0082         target(200+round(rand*PRANGE),:,k)=1;
0083         target(:,:,k)=filter(B,1,target(:,:,k));
0084     end
0085 end
0086 
0087 SNR=0.001;
0088 mixmatrix=randn(1,nchans);
0089 [dummy,maxidx]=max(abs(mixmatrix));
0090 target=nt_mmat(target,mixmatrix);
0091 meg=nt_normcol(meg)+SNR*target;
0092 
0093 
0094 y=nt_normcol(nt_pca(nt_normcol(meg)));
0095 NKEEP=28;
0096 yy=zeros(size(meg,1),NKEEP*(NKEEP+1),size(meg,3));
0097 ii=1;
0098 for k=1:NKEEP;
0099     for j=1:k;
0100         yy(:,ii,:)=meg(:,k,:).*meg(:,j,:);
0101         ii=ii+1;
0102     end
0103 end
0104 
0105 yy=nt_demean2(yy);
0106 SMOOTH=100;
0107 yy=filter(ones(1,SMOOTH),1,yy);
0108 yy=yy(SMOOTH:end,:,:);
0109 
0110 [todss,pwr0,pwr1]=nt_dss1(yy,[],[],0);
0111 z=nt_mmat(yy,todss);
0112 
0113 
0114 if 0
0115     c0=nt_cov(meg(SMOOTH:end,:,:));
0116     c1=nt_cov(meg(SMOOTH:end,:,:).*repmat(z(:,1,:),[1,nchans,1]));
0117 else
0118     %c0=nt_cov([y1,y2]);
0119     c0=nt_cov(meg(SMOOTH:end,:,:).*repmat(max(0,z(:,1,:)),[1,nchans,1]));
0120     c1=nt_cov(meg(SMOOTH:end,:,:).*repmat(min(0,z(:,1,:)),[1,nchans,1]));
0121 end
0122    
0123 [todss2,pwr0,pwr1]=nt_dss0(c0,c1);
0124 z2=nt_mmat(meg,todss2);
0125 
0126 t=(0:size(target,1)-1)/sr;
0127 subplot 234; 
0128 plot(t, squeeze(target(:,1,1:10))/max(max(target(:,1,:))), 'k'); title('target (10 trials)', 'fontsize', 14); set(gca,'ytick',[]); 
0129 xlim([t(1) t(end)]); ylim([-1.2 1.2]);
0130 subplot 235; 
0131 plot(t, squeeze(meg(:,maxidx,:)));  title('mixture', 'fontsize', 14); xlabel('s', 'fontsize',14); set(gca,'ytick',[]); xlim([t(1) t(end)]);
0132 subplot 236; 
0133 plot(t, -squeeze(z2(:,end,1:10))/max(max(abs(z2(:,end,:)))), 'k');  title('recovered', 'fontsize', 14); set(gca,'ytick',[]); 
0134 xlim([t(1) t(end)]); ylim([-1.2 1.2]);
0135 
0136 
0137 figure(2); clf
0138 plot(D);
0139 title('quality of fit to a square');
0140 xlabel('component'); ylabel('score');
0141 
0142 
0143