function mseqq=mseqgen(q,cb)
% Modified M-Sequence of p^n Event Types and a blank
%
% Modified from a m-sequence of p^n Event Types by coding pairs of numbers to
% represent one event class.
%      Think [0 1]=1, [0 -1]=2
%      To keep pairs independent of original cb, m-sequence is shifted by
%      cb.
%
% [mseq2] = mmseq2(powerVal)
%
% q        - the number of event types that are required
%              Condition: Must be p^n for some integer n, and prime, p.
% cb       - highest order of counterbalancing required
%              Condition: Must be less than or equal to powerVal/2
%
% Our base mseq function limits us to p= 2, 3, or 5, as well as what powers
% of these primes we can use.
%
%
% Contact:
% Wesley Kerr
% wesleytk@sas.upenn.edu
%
%
% mseq function inside this is copyrighted with the following formation:
% (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.
 
%% Checking Conditions
 
f=factor(q);
n=length(f);
p=f(1,1);
for i=1:n
    for j=1:n
        if f(1,i)~=f(1,j)
            error 'q must be a power of a single prime'
        end
    end
end
 
if p>5
    error 'Limited by mseq function. Edit mseq function to run this q value.'
end
 
%% Calculating Sequence
 
powerVal=n*(cb+1);
premseq =mseq(p,powerVal);
nlength =length(premseq);
nlength1=nlength-cb;
counter1=1;
mseq2=zeros(nlength1,1);
vals=zeros(1,n);
 
%% Brute Force If-Loopy Way
 
for i=1:nlength1
    shifty= mod(i+3,nlength1)+1;
    if premseq(i)==0
        if premseq(shifty)==0
            mseq2(counter1)=0;
        elseif premseq(shifty)==1
            mseq2(counter1)=1;
        elseif premseq(shifty)==-1
            mseq2(counter1)=2;
        end
    elseif premseq(i)==1
        if premseq(shifty)==0
            mseq2(counter1)=3;
        elseif premseq(shifty)==1
            mseq2(counter1)=4;
        elseif premseq(shifty)==-1
            mseq2(counter1)=-1;
        end
    elseif premseq(i)==-1
        if premseq(shifty)==0
            mseq2(counter1)=-2;
        elseif premseq(shifty)==1
            mseq2(counter1)=-3;
        elseif premseq(shifty)==-1
            mseq2(counter1)=-4;
        end
    end
    counter1=counter1+1;
end
 
% mseq2(727,1)=4;
 
%% New Way -- If you want to create n>2 think about this, but it doesn't
%% work
 
% cbmx=zeros(1,n);
% 
% for i=2:n
%     cbmx(1,i)=mod((i-1)*cb,nlength);
% end
% 
% for i=0:nlength
%     for j=1:n
%         loc(1,j)=mod(i+(j-1)*cb,nlength)+1
%         vals(1,j)=premseq(mod(i+(j-1)*cb,nlength)+1,1);
%         vals=mod(vals,p);
%     end
%     dec=sprintf('%d',vals);
%     mseq2(counter1)=base2dec(dec,p);
%     counter1=counter1+1;
% end
 
 
mseqq=mseq2;