test_nt_quad2square

0001 reset(RandStream.getGlobalStream); % to get reproducible signals
0002 
0003 %{
0004 The data are 10-channel data organized into trials.
0005 The stimulus is a single source consisting of a modulated sinusoid with
0006 random phase. It is mixed into the 10 channels via a 1x10 random mixing matrix.
0007 The noise is produced by 9 independent sources mixed into the 10 channels via a
0008 9x10 random mixing matrix.
0009 %}
0010 
0011 sr=1000;
0012 nsamples=1000;
0013 ntrials=100;
0014 nchans=10;
0015 
0016 % sinusoidal pulse with random phase
0017 CF1=4;
0018 target=sin(2*pi*0.5*(1:nsamples)'/sr).^2;
0019 target=repmat(target,[1,1,ntrials]) .* ...
0020     sin( 2*pi*(CF1*repmat((1:nsamples)',[1,1,ntrials])/sr + ...
0021     repmat(rand(1,1,ntrials),[nsamples,1,1]))); % phase
0022 target=nt_normcol(target);
0023 
0024 % noise
0025 NNOISE=9;
0026 noise=randn(nsamples,NNOISE,ntrials);
0027 noise=nt_mmat(noise,randn(NNOISE,nchans));
0028 noise=nt_normcol(noise);
0029 
0030 % data
0031 SNR=0.0001;
0032 x=noise+ SNR * nt_mmat(target,randn(1,nchans)) ;
0033 x=nt_demean(x);
0034 x=nt_normcol(x);
0035 
0036 
0037 %{
0038 We append a DC channel (non-zero constant value) and then we form all
0039 cross-products of channels two-by-two.  Appending a DC channel implies
0040 that the set of cross products also includes the original data channels, 
0041 as well as the DC channel.
0042 %}
0043 
0044 % append DC channel
0045 x=[x, ones(nsamples,1,ntrials)*mean(abs(x(:)))];
0046 
0047 % not sure why this makes things better:
0048 if 1
0049     THRESH=0;%10^-12;
0050     [topcs,pwr]=nt_pca0(x,[],[],THRESH); 
0051     x=nt_mmat(x,topcs);
0052     x=nt_normcol(x);
0053 end
0054 
0055 % form cross products
0056 xx=nt_xprod(x);
0057 
0058 % DSS finds the most reproducible quadratic form
0059 keep1=[]; keep2=0;%10^-12;
0060 [toquads,pwr0,pwr1]=nt_dss1(xx,[],keep1,keep2);
0061 
0062 % illustrate
0063 figure(1); clf
0064 plot(pwr1./pwr0,'.-')
0065 z=nt_mmat(xx,toquads);
0066 figure(2); clf
0067 subplot 211; nt_bsplot(z(:,1,:));  title('DC')
0068 subplot 212; nt_bsplot(z(:,2,:)); title('most reproducible quadratic form');
0069 
0070 
0071 %{
0072 We now find the linear component (weighted sum of channels) with square
0073 closest to our optimal quadratic form.
0074 %}
0075 
0076 [tosquares,D]=nt_quad2square(toquads(:,2));
0077 z=nt_mmat(x,tosquares);
0078 
0079 figure(3); clf
0080 subplot 211; nt_bsplot(z(:,1,:).^2); 
0081 ylabel('power'); title('closest square');
0082 subplot 212; nt_bsplot(z(:,2,:).^2); 
0083 ylabel('power'); title('second closest'); 
0084 figure(4); clf
0085 plot(abs(D(2:end)))
0086 
0087 %{
0088 It may necessary to sort the components to ensure that the DC component
0089 does not come first.
0090 %}
0091