module OW_mseq
!		  Maximum length sequence assuming 2,3,5 distinct values
!
!       [ms]=MSEQ(baseVal,powerVal[,shift,whichSeq])
!
!       OUTPUT:
!       ms = generated maximum length sequence, of length basisVal^powerVal-1
!
!       INPUT:
!		  baseVal  -nuber of sequence levels (2,3, or 5 allowed)
!		  powerVal -power, so that sequence length is baseVal^powerVal-1
!		  shift    -cyclical shift of the sequence
!		  whichSeq -sequence istantiation to use 
!		  (numer of sequences varies with powreVal - see the code)

! (c) Giedrius T. Buracas, SNL-B, Salk Institute
! Register values are taken from: WDT Davies, System Identification
! for self-adaptive control. Wiley-Interscience, 1970
! When using mseq code for design of FMRI experiments, please, cite:
! G.T.Buracas & G.M.Boynton (2002) Efficient Design of Event-Related fMRI 
! Experiments Using M-sequences. NeuroImage, 16, 801-813.
 

subroutine owmseq(baseVal,powerVal,shift,whichSeq)

optional, integer				:: shift,whichSeq
integer 					:: baseVal,powerVal
integer					:: bitNum
integer, allocatable			:: register(:)
real					:: ms(:)

if (not(present(whichSeq))) whichSeq=1
if (not(present(shift)))    shift=1

bitNum=baseVal^powerVal-1
allocate(register(powerVal)); register=1

if (powerVal .le. 30) then
	if (baseVal==2) then
		select case (powerVal)
			case (2); tap(1).No=[1,2];
			case (3); tap(1).No=[1,3];
			 	tap(2).No=[2,3];
			case (4); tap(1).No=[1,4];
			 	tap(2).No=[3,4];
			case (5); tap(1).No=[2,5];
				tap(2).No=[3,5];
				tap(3).No=[1,2,3,5];
				tap(4).No=[2,3,4,5];
				tap(5).No=[1,2,4,5];
				tap(6).No=[1,3,4,5];
			case (6); tap(1).No=[1,6];
				tap(2).No=[5,6];
				tap(3).No=[1,2,5,6];
				tap(4).No=[1,4,5,6];
				tap(5).No=[1,3,4,6];
				tap(6).No=[2,3,5,6];
			case (7); tap(1).No=[1,7];
				tap(2).No=[6,7];
				tap(3).No=[3,7];
				tap(4).No=[4,7];
				tap(5).No=[1,2,3,7];
				tap(6).No=[4,5,6,7];
				tap(7).No=[1,2,5,7];
				tap(8).No=[2,5,6,7];
				tap(9).No=[2,3,4,7];
				tap(10).No=[3,4,5,7];
				tap(11).No=[1,3,5,7];
				tap(12).No=[2,4,6,7];
				tap(13).No=[1,3,6,7];
				tap(14).No=[1,4,6,7];
				tap(15).No=[2,3,4,5,6,7];
				tap(16).No=[1,2,3,4,5,7];
				tap(17).No=[1,2,4,5,6,7];
				tap(18).No=[1,2,3,5,6,7];
			case (8); tap(1).No=[1,2,7,8];
				tap(2).No=[1,6,7,8];
				tap(3).No=[1,3,5,8];
				tap(4).No=[3,5,7,8];
				tap(5).No=[2,3,4,8];
				tap(6).No=[4,5,6,8];
				tap(7).No=[2,3,5,8];
				tap(8).No=[3,5,6,8];
				tap(9).No=[2,3,6,8];
				tap(10).No=[2,5,6,8];
				tap(11).No=[2,3,7,8];
				tap(12).No=[1,5,6,8];
				tap(13).No=[1,2,3,4,6,8];
				tap(14).No=[2,4,5,6,7,8];
				tap(15).No=[1,2,3,6,7,8];
				tap(16).No=[1,2,5,6,7,8];
			case (9); tap(1).No=[4,9];
				tap(2).No=[5,9];
				tap(3).No=[3,4,6,9];
				tap(4).No=[3,5,6,9];
				tap(5).No=[4,5,8,9];
				tap(6).No=[1,4,5,9];
				tap(7).No=[1,4,8,9];
				tap(8).No=[1,5,8,9];
				tap(9).No=[2,3,5,9];
				tap(10).No=[4,6,7,9];
				tap(11).No=[5,6,8,9];
				tap(12).No=[1,3,4,9];
				tap(13).No=[2,7,8,9];
				tap(14).No=[1,2,7,9];
				tap(15).No=[2,4,7,9];
				tap(16).No=[2,5,7,9];
				tap(17).No=[2,4,8,9];
				tap(18).No=[1,5,7,9];
				tap(19).No=[1,2,4,5,6,9];
				tap(20).No=[3,4,5,7,8,9];
				tap(21).No=[1,3,4,6,7,9];
				tap(22).No=[2,3,5,6,8,9];
				tap(23).No=[3,5,6,7,8,9];
				tap(24).No=[1,2,3,4,6,9];
				tap(25).No=[1,5,6,7,8,9];
				tap(26).No=[1,2,3,4,8,9];
				tap(27).No=[1,2,3,7,8,9];
				tap(28).No=[1,2,6,7,8,9];
				tap(29).No=[1,3,5,6,8,9];
				tap(30).No=[1,3,4,6,8,9];
				tap(31).No=[1,2,3,5,6,9];
				tap(32).No=[3,4,6,7,8,9];
				tap(33).No=[2,3,6,7,8,9];
				tap(34).No=[1,2,3,6,7,9];
				tap(35).No=[1,4,5,6,8,9];
				tap(36).No=[1,3,4,5,8,9];
				tap(37).No=[1,3,6,7,8,9];
				tap(38).No=[1,2,3,6,8,9];
				tap(39).No=[2,3,4,5,6,9];
				tap(40).No=[3,4,5,6,7,9];
				tap(41).No=[2,4,6,7,8,9];
				tap(42).No=[1,2,3,5,7,9];
				tap(43).No=[2,3,4,5,7,9];
				tap(44).No=[2,4,5,6,7,9];
				tap(45).No=[1,2,4,5,7,9];
				tap(46).No=[2,4,5,6,7,9];
				tap(47).No=[1,3,4,5,6,7,8,9];
				tap(48).No=[1,2,3,4,5,6,8,9];
			case (10);tap(1).No=[3,10];
				tap(2).No=[7,10];
				tap(3).No=[2,3,8,10];
				tap(4).No=[2,7,8,10];
				tap(5).No=[1,3,4,10];
				tap(6).No=[6,7,9,10];
				tap(7).No=[1,5,8,10];
				tap(8).No=[2,5,9,10];
				tap(9).No=[4,5,8,10];
				tap(10).No=[2,5,6,10];
				tap(11).No=[1,4,9,10];
				tap(12).No=[1,6,9,10];
				tap(13).No=[3,4,8,10];
				tap(14).No=[2,6,7,10];
				tap(15).No=[2,3,5,10];
				tap(16).No=[5,7,8,10];
				tap(17).No=[1,2,5,10];
				tap(18).No=[5,8,9,10];
				tap(19).No=[2,4,9,10];
				tap(20).No=[1,6,8,10];
				tap(21).No=[3,7,9,10];
				tap(22).No=[1,3,7,10];
				tap(23).No=[1,2,3,5,6,10];
				tap(24).No=[4,5,7,8,9,10];
				tap(25).No=[2,3,6,8,9,10];
				tap(26).No=[1,2,4,7,8,10];
				tap(27).No=[1,5,6,8,9,10];
				tap(28).No=[1,2,4,5,9,10];
				tap(29).No=[2,5,6,7,8,10];
				tap(30).No=[2,3,4,5,8,10];
				tap(31).No=[2,4,6,8,9,10];
				tap(32).No=[1,2,4,6,8,10];
				tap(33).No=[1,2,3,7,8,10];
				tap(34).No=[2,3,7,8,9,10];
				tap(35).No=[3,4,5,8,9,10];
				tap(36).No=[1,2,5,6,7,10];
				tap(37).No=[1,4,6,7,9,10];
				tap(38).No=[1,3,4,6,9,10];
				tap(39).No=[1,2,6,8,9,10];
				tap(40).No=[1,2,4,8,9,10];
				tap(41).No=[1,4,7,8,9,10];
				tap(42).No=[1,2,3,6,9,10];
				tap(43).No=[1,2,6,7,8,10];
				tap(44).No=[2,3,4,8,9,10];
				tap(45).No=[1,2,4,6,7,10];
				tap(46).No=[3,4,6,8,9,10];
				tap(47).No=[2,4,5,7,9,10];
				tap(48).No=[1,3,5,6,8,10];
				tap(49).No=[3,4,5,6,9,10];
				tap(50).No=[1,4,5,6,7,10];
				tap(51).No=[1,3,4,5,6,7,8,10];
				tap(52).No=[2,3,4,5,6,7,9,10];
				tap(53).No=[3,4,5,6,7,8,9,10];
				tap(54).No=[1,2,3,4,5,6,7,10];
				tap(55).No=[1,2,3,4,5,6,9,10];
				tap(56).No=[1,4,5,6,7,8,9,10];
				tap(57).No=[2,3,4,5,6,8,9,10];
				tap(58).No=[1,2,4,5,6,7,8,10];
				tap(59).No=[1,2,3,4,6,7,9,10];
				tap(60).No=[1,3,4,6,7,8,9,10];
			case (11); tap(1).No=[9,11];
			case (12); tap(1).No=[6,8,11,12];
			case (13); tap(1).No=[9,10,12,13];
			case (14); tap(1).No=[4,8,13,14];
			case (15); tap(1).No=[14,15];
			case (16); tap(1).No=[4,13,15,16];
			case (17); tap(1).No=[14,17];
			case (18); tap(1).No=[11,18];
			case (19); tap(1).No=[14,17,18,19];
			case (20); tap(1).No=[17,20];
			case (21); tap(1).No=[19,21];
			case (22); tap(1).No=[21,22];
			case (23); tap(1).No=[18,23];
			case (24); tap(1).No=[17,22,23,24];
			case (25); tap(1).No=[22,25];
			case (26); tap(1).No=[20,24,25,26];
			case (27); tap(1).No=[22,25,26,27];
			case (28); tap(1).No=[25,28];
			case (29); tap(1).No=[27,29];
			case (30); tap(1).No=[7,28,29,30];
		end select
	else if (baseVal==3) then
		select case (powerVal)
			case (2);	   tap(1).No=[2,1];
				   tap(2).No=[1,1];
			case (3);	   tap(1).No=[0,1,2];
				   tap(2).No=[1,0,2];
				   tap(3).No=[1,2,2];
				   tap(4).No=[2,1,2];
			case (4);	   tap(1).No=[0,0,2,1];
				   tap(2).No=[0,0,1,1];
				   tap(3).No=[2,0,0,1];
				   tap(4).No=[2,2,1,1];
				   tap(5).No=[2,1,1,1];
				   tap(6).No=[1,0,0,1];
				   tap(7).No=[1,2,2,1];
				   tap(8).No=[1,1,2,1];
			case (5);	   tap(1).No=[0,0,0,1,2]; 
				   tap(2).No=[0,0,0,1,2];
				   tap(3).No=[0,0,1,2,2];
				   tap(4).No=[0,2,1,0,2];
				   tap(5).No=[0,2,1,1,2];
				   tap(6).No=[0,1,2,0,2];
				   tap(7).No=[0,1,1,2,2];
				   tap(8).No=[2,0,0,1,2];
				   tap(9).No=[2,0,2,0,2];
				   tap(10).No=[2,0,2,2,2];
				   tap(11).No=[2,2,0,2,2];
				   tap(12).No=[2,2,2,1,2];
				   tap(13).No=[2,2,1,2,2];
				   tap(14).No=[2,1,2,2,2];
				   tap(15).No=[2,1,1,0,2];
				   tap(16).No=[1,0,0,0,2];
				   tap(17).No=[1,0,0,2,2];
				   tap(18).No=[1,0,1,1,2];
				   tap(19).No=[1,2,2,2,2];
				   tap(20).No=[1,1,0,1,2];
				   tap(21).No=[1,1,2,0,2];
			case (6);	   tap(1).No=[0,0,0,0,2,1];
				   tap(2).No=[0,0,0,0,1,1];
				   tap(3).No=[0,0,2,0,2,1];
				   tap(4).No=[0,0,1,0,1,1];
				   tap(5).No=[0,2,0,1,2,1];
				   tap(6).No=[0,2,0,1,1,1];
				   tap(7).No=[0,2,2,0,1,1];
				   tap(8).No=[0,2,2,2,1,1];
				   tap(9).No=[2,1,1,1,0,1];
				   tap(10).No=[1,0,0,0,0,1];
				   tap(11).No=[1,0,2,1,0,1];
				   tap(12).No=[1,0,1,0,0,1];
				   tap(13).No=[1,0,1,2,1,1];
				   tap(14).No=[1,0,1,1,1,1];
				   tap(15).No=[1,2,0,2,2,1];
				   tap(16).No=[1,2,0,1,0,1];
				   tap(17).No=[1,2,2,1,2,1];
				   tap(18).No=[1,2,1,0,1,1];
				   tap(19).No=[1,2,1,2,1,1];
				   tap(20).No=[1,2,1,1,2,1];
				   tap(21).No=[1,1,2,1,0,1];
				   tap(22).No=[1,1,1,0,1,1];
				   tap(23).No=[1,1,1,2,0,1];
				   tap(24).No=[1,1,1,1,1,1];
			case (7);	   tap(1).No=[0,0,0,0,2,1,2];
				   tap(2).No=[0,0,0,0,1,0,2];
				   tap(3).No=[0,0,0,2,0,2,2];
				   tap(4).No=[0,0,0,2,2,2,2];
				   tap(5).No=[0,0,0,2,1,0,2];
				   tap(6).No=[0,0,0,1,1,2,2];
				   tap(7).No=[0,0,0,1,1,1,2];
				   tap(8).No=[0,0,2,2,2,0,2];
				   tap(9).No=[0,0,2,2,1,2,2];
				   tap(10).No=[0,0,2,1,0,0,2];
				   tap(11).No=[0,0,2,1,2,2,2];
				   tap(12).No=[0,0,1,0,2,1,2];
				   tap(13).No=[0,0,1,0,1,1,2];
				   tap(14).No=[0,0,1,1,0,1,2];
				   tap(15).No=[0,0,1,1,2,0,2];
				   tap(16).No=[0,2,0,0,0,2,2];
				   tap(17).No=[0,2,0,0,1,0,2];
				   tap(18).No=[0,2,0,0,1,1,2];
				   tap(19).No=[0,2,0,2,2,0,2];
				   tap(20).No=[0,2,0,2,1,2,2];
				   tap(21).No=[0,2,0,1,1,0,2];
				   tap(22).No=[0,2,2,0,2,0,2];
				   tap(23).No=[0,2,2,0,1,2,2];
				   tap(24).No=[0,2,2,2,2,1,2];
				   tap(25).No=[0,2,2,2,1,0,2];
				   tap(26).No=[0,2,2,1,0,1,2];
				   tap(27).No=[0,2,2,1,2,2,2];
		end select   
elseif baseVal==5,
switch powerVal,
case 2, tap(1).No=[4,3];
		   tap(2).No=[3,2];
		   tap(3).No=[2,2];
		   tap(4).No=[1,3];
case 3, tap(1).No=[0,2,3];
		   tap(2).No=[4,1,2];
		   tap(3).No=[3,0,2];
		   tap(4).No=[3,4,2];
		   tap(5).No=[3,3,3];
		   tap(6).No=[3,3,2];
		   tap(7).No=[3,1,3];
		   tap(8).No=[2,0,3];
		   tap(9).No=[2,4,3];
		   tap(10).No=[2,3,3];
		   tap(11).No=[2,3,2];
		   tap(12).No=[2,1,2];
		   tap(13).No=[1,0,2];
		   tap(14).No=[1,4,3];
		   tap(15).No=[1,1,3];
case 4, tap(1).No=[0,4,3,3];
		   tap(2).No=[0,4,3,2];
		   tap(3).No=[0,4,2,3];
		   tap(4).No=[0,4,2,2];
		   tap(5).No=[0,1,4,3];
		   tap(6).No=[0,1,4,2];
		   tap(7).No=[0,1,1,3];
		   tap(8).No=[0,1,1,2];
		   tap(9).No=[4,0,4,2];
		   tap(10).No=[4,0,3,2];
		   tap(11).No=[4,0,2,3];
		   tap(12).No=[4,0,1,3];
		   tap(13).No=[4,4,4,2];
		   tap(14).No=[4,3,0,3];
		   tap(15).No=[4,3,4,3];
		   tap(16).No=[4,2,0,2];
		   tap(17).No=[4,2,1,3];
		   tap(18).No=[4,1,1,2];
		   tap(19).No=[3,0,4,2];
		   tap(20).No=[3,0,3,3];
		   tap(21).No=[3,0,2,2];
		   tap(22).No=[3,0,1,3];
		   tap(23).No=[3,4,3,2];
		   tap(24).No=[3,3,0,2];
		   tap(25).No=[3,3,3,3];
		   tap(26).No=[3,2,0,3];
		   tap(27).No=[3,2,2,3];
		   tap(28).No=[3,1,2,2];
		   tap(29).No=[2,0,4,3];
		   tap(30).No=[2,0,3,2];
		   tap(31).No=[2,0,2,3];
		   tap(32).No=[2,0,1,2];
		   tap(33).No=[2,4,2,2];
		   tap(34).No=[2,3,0,2];
		   tap(35).No=[2,3,2,3];
		   tap(36).No=[2,2,0,3];
		   tap(37).No=[2,2,3,3];
		   tap(38).No=[2,1,3,2];
		   tap(39).No=[1,0,4,3];
		   tap(40).No=[1,0,3,3];
		   tap(41).No=[1,0,2,2];
		   tap(42).No=[1,0,1,2];
		   tap(43).No=[1,4,1,2];
		   tap(44).No=[1,3,0,3];
		   tap(45).No=[1,3,1,3];
		   tap(46).No=[1,2,0,2];
		   tap(47).No=[1,2,4,3];
		   tap(48).No=[1,1,4,2];
end if
else
call seperr ("M-sequence is not defined for powerVal > 30"))
end if

ms=zeros(bitNum,1);
if isempty(whichSeq), whichSeq=ceil(rand(1)*length(tap)); 
else, 
   if whichSeq>length(tap) | whichSeq<1
      disp(sprintf(' wrapping arround!'));
      whichSeq=rem(whichSeq,length(tap))+1;
   end;
end;

weights=zeros(1,powerVal);
if baseVal==2,
	weights(tap(whichSeq).No)=1;
elseif baseVal>2,
   weights=tap(whichSeq).No;
end;
  
!weights
  
for i=1:bitNum
   ! calculating next digit with modulo powerVal arithmetic
   !   register, (tap(1).No)

  !ms(i)=mod(sum(register(tap(whichSeq).No)),baseVal);
  ms(i)=mod(weights*register+baseVal,baseVal);
  ! updating the register
  
  register=[ms(i);register(1:powerVal-1)];
end

ms=ms(:);
if ~isempty(shift),
   shift=rem(shift, length(ms));
   ms=[ms(shift+1:end); ms(1:shift)];
end;

if baseVal==2,     
  ms=ms*2-1;
elseif baseVal==3, 
  ms(ms==2)=-1;
elseif baseVal==5, 
  ms(ms==4)=-1;
  ms(ms==3)=-2;
else
  error('wrong baseVal!');
end;


