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(:) integer, dimension (1:179) :: index=(/2, & 4, & 6, & 8, & 10, & 12, & 14, & 18, & 22, & 26, & 30, & 32, & 34, & 38, & 42, & 46, & 50, & 52, & 54, & 56, & 58, & 62, & 66, & 70, & 74, & 78, & 82, & 86, & 90, & 94, & 98, & 104, & 110, & 116, & 120, & 124, & 128, & 132, & 136, & 140, & 144, & 148, & 152, & 156, & 160, & 164, & 168, & 174, & 180, & 186, & 190, & 192, & 194, & 198, & 202, & 206, & 210, & 214, & 218, & 222, & 226, & 230, & 234, & 238, & 242, & 246, & 250, & 254, & 258, & 264, & 270, & 276, & 282, & 288, & 294, & 300, & 306, & 312, & 318, & 324, & 330, & 336, & 342, & 348, & 354, & 360, & 366, & 372, & 378, & 384, & 390, & 396, & 402, & 408, & 414, & 420, & 426, & 434, & 442, & 444, & 446, & 450, & 454, & 458, & 462, & 466, & 470, & 474, & 478, & 482, & 486, & 490, & 494, & 498, & 502, & 506, & 510, & 514, & 518, & 522, & 526, & 532, & 538, & 544, & 550, & 556, & 562, & 568, & 574, & 580, & 586, & 592, & 598, & 604, & 610, & 616, & 622, & 628, & 634, & 640, & 646, & 652, & 658, & 664, & 670, & 676, & 682, & 688, & 694, & 702, & 710, & 718, & 726, & 734, & 742, & 750, & 758, & 766, & 774, & 776, & 780, & 784, & 788, & 790, & 794, & 796, & 798, & 802, & 804, & 806, & 808, & 810, & 814, & 816, & 820, & 824, & 826, & 828, & 832 /) integer, dimension (1:179) :: seq = (/ 1,2, & 1,3, & 2,3, & 1,4, & 3,4, & 2,5, & 3,5, & 1,2,3,5, & 2,3,4,5, & 1,2,4,5, & 1,3,4,5, & 1,6, & 5,6, & 1,2,5,6, & 1,4,5,6, & 1,3,4,6, & 2,3,5,6, & 1,7, & 6,7, & 3,7, & 4,7, & 1,2,3,7, & 4,5,6,7, & 1,2,5,7, & 2,5,6,7, & 2,3,4,7, & 3,4,5,7, & 1,3,5,7, & 2,4,6,7, & 1,3,6,7, & 1,4,6,7, & 2,3,4,5,6,7, & 1,2,3,4,5,7, & 1,2,4,5,6,7, & 1,2,3,5,6,7, & 1,2,7,8, & 1,6,7,8, & 1,3,5,8, & 3,5,7,8, & 2,3,4,8, & 4,5,6,8, & 2,3,5,8, & 3,5,6,8, & 2,3,6,8, & 2,5,6,8, & 2,3,7,8, & 1,5,6,8, & 1,2,3,4,6,8, & 2,4,5,6,7,8, & 1,2,3,6,7,8, & 1,2,5,6,7,8, & 4,9, & 5,9, & 3,4,6,9, & 3,5,6,9, & 4,5,8,9, & 1,4,5,9, & 1,4,8,9, & 1,5,8,9, & 2,3,5,9, & 4,6,7,9, & 5,6,8,9, & 1,3,4,9, & 2,7,8,9, & 1,2,7,9, & 2,4,7,9, & 2,5,7,9, & 2,4,8,9, & 1,5,7,9, & 1,2,4,5,6,9, & 3,4,5,7,8,9, & 1,3,4,6,7,9, & 2,3,5,6,8,9, & 3,5,6,7,8,9, & 1,2,3,4,6,9, & 1,5,6,7,8,9, & 1,2,3,4,8,9, & 1,2,3,7,8,9, & 1,2,6,7,8,9, & 1,3,5,6,8,9, & 1,3,4,6,8,9, & 1,2,3,5,6,9, & 3,4,6,7,8,9, & 2,3,6,7,8,9, & 1,2,3,6,7,9, & 1,4,5,6,8,9, & 1,3,4,5,8,9, & 1,3,6,7,8,9, & 1,2,3,6,8,9, & 2,3,4,5,6,9, & 3,4,5,6,7,9, & 2,4,6,7,8,9, & 1,2,3,5,7,9, & 2,3,4,5,7,9, & 2,4,5,6,7,9, & 1,2,4,5,7,9, & 2,4,5,6,7,9, & 1,3,4,5,6,7,8,9, & 1,2,3,4,5,6,8,9, & 3,10, & 7,10, & 2,3,8,10, & 2,7,8,10, & 1,3,4,10, & 6,7,9,10, & 1,5,8,10, & 2,5,9,10, & 4,5,8,10, & 2,5,6,10, & 1,4,9,10, & 1,6,9,10, & 3,4,8,10, & 2,6,7,10, & 2,3,5,10, & 5,7,8,10, & 1,2,5,10, & 5,8,9,10, & 2,4,9,10, & 1,6,8,10, & 3,7,9,10, & 1,3,7,10, & 1,2,3,5,6,10, & 4,5,7,8,9,10, & 2,3,6,8,9,10, & 1,2,4,7,8,10, & 1,5,6,8,9,10, & 1,2,4,5,9,10, & 2,5,6,7,8,10, & 2,3,4,5,8,10, & 2,4,6,8,9,10, & 1,2,4,6,8,10, & 1,2,3,7,8,10, & 2,3,7,8,9,10, & 3,4,5,8,9,10, & 1,2,5,6,7,10, & 1,4,6,7,9,10, & 1,3,4,6,9,10, & 1,2,6,8,9,10, & 1,2,4,8,9,10, & 1,4,7,8,9,10, & 1,2,3,6,9,10, & 1,2,6,7,8,10, & 2,3,4,8,9,10, & 1,2,4,6,7,10, & 3,4,6,8,9,10, & 2,4,5,7,9,10, & 1,3,5,6,8,10, & 3,4,5,6,9,10, & 1,4,5,6,7,10, & 1,3,4,5,6,7,8,10, & 2,3,4,5,6,7,9,10, & 3,4,5,6,7,8,9,10, & 1,2,3,4,5,6,7,10, & 1,2,3,4,5,6,9,10, & 1,4,5,6,7,8,9,10, & 2,3,4,5,6,8,9,10, & 1,2,4,5,6,7,8,10, & 1,2,3,4,6,7,9,10, & 1,3,4,6,7,8,9,10, & 9,11, & 6,8,11,12, & 9,10,12,13, & 4,8,13,14, & 14,15, & 4,13,15,16, & 14,17, & 11,18, & 14,17,18,19, & 17,20, & 19,21, & 21,22, & 18,23, & 17,22,23,24, & 22,25, & 20,24,25,26, & 22,25,26,27, & 25,28, & 27,29, & 7,28,29,30/) if (not(present(whichSeq))) whichSeq=1 if (not(present(shift))) shift=1 baseVal = 2 bitNum=baseVal^powerVal-1 allocate(register(powerVal)); register=1 lengtap=1 if (powerVal .le. 30) then select case (powerVal) case (3) lengtap=2 case (4) lengtap=2 case (5) lengtap=6 case (6) lengtap=6 case (7) lengtap=18 case (8) lengtap=16 case (9) lengtap=48 case (10) lengtap=60 end select allocate(ms(bitNum)); ms=0. if (not(present).whichSeq) then call owexec_time (rnd) seed(1)=int(rnd/10.) call owexec_time (rnd) seed(2)=int(rnd+2000.) call random_seed(put=seed) call random_number(rnd) whichSeq=ceiling(rnd*lengtap) else if (whichSeq > lengtap .or. whichSeq<1) then write(0,*) ' WARNING: wrapping arround!' whichSeq=mod(whichSeq,lengtap)+1; end if; end if; allocate(weights(powerVal)); weights=0. weights(tap(whichSeq).No)=1; for i=1:bitNum ! calculating next digit with modulo powerVal arithmetic ! register, (tap(1).No) ms(i)=mod(weights*register+baseVal,baseVal); ! updating the register register=[ms(i);register(1:powerVal-1)]; end do ms=ms(:); if ~isempty(shift), shift=mod(shift, length(ms)); ms=[ms(shift+1:end); ms(1:shift)]; end if if baseVal==2, ms=ms*2-1; else if baseVal==3, ms(ms==2)=-1; else if baseVal==5, ms(ms==4)=-1; ms(ms==3)=-2; else call seperr("wrong baseVal!"); end if else call seperr ("M-sequence is not defined for powerVal > 30")) end if