function [data]=simuldata(N,P,gamma,myswitch) %generates one example for the causality challenge % % Input: % N : number of time points % P : order of AR-system % gamma: parameter controlling the relative strength of noise % and signal, gamma=0 means only signal, gamma=1 means only noise % myswitch: if positive direction is from channel 1 to channel 2; % if negative direction is from channel 2 to channel 1; % % Output: % data Nx2 matrix of data % % One example of the challenge was generated with the command: % N=6000; P=10; gamma=rand;myswitch=randn; % data=simuldata(N,P,gamma,myswitch) % ssig=1-gamma; snoise=gamma; nrun=1; M=2;M1=2; M2=3; ckount=0; while ckount< nrun Arsig=[]; for k=1:P aloc=randn(M)/4; if myswitch<0; aloc(2,1)=0; else aloc(1,2)=0; end Arsig=[Arsig,aloc]; end E=eye(M*P);AA=[Arsig;E(1:end-M,:)];lambda=eig(AA);lambdamax=max(abs(lambda)); Arnoise=[]; for k=1:P;aloc=diag(diag(randn(M2)/4)); Arnoise=[Arnoise,aloc];end E=eye(M2*P);AAnoise=[Arnoise;E(1:end-M2,:)]; lambda=eig(AAnoise); lambdamaxnoise=max(abs(lambda)); %subplot(3,3,i);plot(imag(ps)); if lambdamax<1 & lambdamaxnoise<1; ckount=ckount+1; x=rand(M1,N)-.5; %x=randn(M1,N); data= mymvfilter(Arsig,x); datasig=data'; siglevel=norm(datasig,'fro'); x=rand(M2,N)-.5; %x=randn(M2,N); data= mymvfilter(Arnoise,x); datanoise=data'; B=randn(M2,M1);datanoise=datanoise*B; noiselevel=norm(datanoise,'fro'); data=ssig*datasig/siglevel+snoise*datanoise/noiselevel; end end return; function y=mymvfilter(Ar,x); [nchan,norder]=size(Ar); norder=norder/nchan; Ar_reshape=zeros(nchan,nchan,norder); for i=1:norder; Ar_reshape(:,:,i)=Ar(:,(i-1)*nchan+1:i*nchan); end [nchan,N]=size(x); y=x; for i=2:N; norderloc=min([i-1,norder]); for j=1:norderloc; y(:,i)=y(:,i)+Ar_reshape(:,:,j)*y(:,i-j); end end return