%   Pgss.m

%   This software is protected by german copyright and international treaties.             
%   Copyright 2004 Markus Junghöfer & Peter Peyk. All Rights Reserved.                     
%                                                                                          
%   THIS SOFTWARE AND DOCUMENTATION IS PROVIDED "AS IS," AND COPYRIGHT HOLDERS MAKE        
%   NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO,    
%   WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR THAT THE USE    
%   OF THE SOFTWARE OR DOCUMENTATION WILL NOT INFRINGE ANY THIRD PARTY PATENTS, COPYRIGHTS,
%   TRADEMARKS OR OTHER RIGHTS. COPYRIGHT HOLDERS WILL NOT BE LIABLE FOR ANY DIRECT,       
%   INDIRECT, SPECIAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF ANY USE OF THE SOFTWARE      
%   OR DOCUMENTATION. 


global xGX

for GroupInd=1:1
    NCond=5
    SPMcfg_path=['D:\ValAro\00_Stats_NoScaling\FixStats\All\SPMcfg.mat']
    SPMcfg_path=['C:\user\Markus\PGSS\All\SPMcfg.mat']
    load(SPMcfg_path)
    if GroupInd==1
        SubIndVec=[1:21]
    elseif GroupInd==2
        FemaleVec=[7 9 10 11 13 14 16 17 18 20 22]
        SubIndVec=[5 6 7 8 10 11 13 14 15 16 18]
    elseif GroupInd==3
        MaleVec =[3 4  5  6  12 15 21 23 24 25] % 2
        SubIndVec=[1 2 3 4 9 12 17 19 20 21]
    end
    NSub=21;
    NSes=2;
    
    Tmp0=find(xX.X(:,5)>1.7);
    Tmp00=[];
    for k=1:length(Tmp0)
        Tmp00=[Tmp00,[Tmp0(k)-2:Tmp0(k)+1]];
    end
    w=xX.X(Tmp00,5);
    MeanMat=zeros(NSub,NSes,NCond,length(w));
    RGVec=xGX.rg;
    SampRate=1/3;
    N=6;
    Nyquist=SampRate./2;
    
    [B,A]=butter(N,2/3,'low');
    FiltRGVec=(filtfilt(B,A,RGVec));
    [B,A]=butter(N,1/10,'high');
    FiltRGVec=(filtfilt(B,A,FiltRGVec))+1;
    %FiltRGVec=RGVec;
    % figure
    % plot(RGVec-mean(RGVec))
    % hold on
    % plot(FiltRGVec,'r')
    % hold off
    for SubInd=1:NSub
        for SesInd=1:NSes
            for i=1:NCond
                iVec=[(SesInd-1).*11+i:22:460];
                Tmp1=find(xX.X(:,iVec(SubInd))>1.7);
                Tmp2=[];
                for k=1:length(Tmp1)
                    Tmp2=[Tmp2,[Tmp1(k)-2:Tmp1(k)+1]];
                end
                MeanMat(SubInd,SesInd,i,:)=FiltRGVec(Tmp2).*w;
            end
        end
    end
    WeightVec=[3 -2 -2 -2 3];
    WeightVec=[1 0 -2 0 1];
    WeightInd=[1 3 5]
%     WeightVec=[0 1 -2 1 0];
%     WeightInd=[2 3 4]
    Weights=WeightVec(WeightInd)
    NCond=length(WeightInd)
    MeanMat=MeanMat(SubIndVec,:,WeightInd,:);
    NSub=length(SubIndVec);
    
    
    a=squeeze(mean(MeanMat,4));
    Fast=squeeze(a(:,1,:));
    Slow=squeeze(a(:,2,:));
    SlowFast=squeeze(mean(a,2));
    %SlowFast=[Slow;Fast];
    %NSub=NSub.*2;
    b=squeeze(mean(a,1));
    c=mean(b);
    
    %     figure
    %     subplot(3,1,1)
    %     plot(b(1,:))
    %     subplot(3,1,2)
    %     plot(b(2,:),'r')
    %     subplot(3,1,3)
    %     plot(c,'g')
    
    figure
    subplot(3,1,1)
    plot(b(1,:)./mean(b(1,:)))
    axis([0 6 0.95 1.05])
    subplot(3,1,2)
    plot(b(2,:)./mean(b(2,:)),'r')
    axis([0 6 0.95 1.05])
    subplot(3,1,3)
    plot(c./mean(c),'g')
    axis([0 6 0.95 1.05])
    

    RepMat=SlowFast;

    %Berechnung von Mean square contrast
    sumrep = sum(RepMat);%summiert ueber Vps   
    L2 = (sum(sumrep.*Weights)).^2;     
    MScont = L2/(NSub*sum(Weights.^2));
    
    %Corrections of data matrix for Grand Mean     
    
    corGM=RepMat-mean(mean(RepMat));
    
    %for row effect    
    rowef=mean(corGM,2);   
    corGMR=zeros(size(RepMat));
    for CondInd=1:NCond
        corGMR(1:NSub,CondInd) = corGM(1:NSub,CondInd)-rowef;
    end  
    %for column effect 
    colef = mean(corGMR);
    corGMRC=zeros(size(RepMat));
    for SubInd =1:NSub
        corGMRC(SubInd,1:NCond)=corGMR(SubInd,1:NCond)-colef;
    end
    
    %sum of squares for condition x subjects
    %this is used as error term for calculation of F. one also could
    %calculate seperate error terms for each trend type. see Rosnow &
    %Rosenthal for details. the approach here is sufficient in most cases.
    
    SScs = sum(sum(corGMRC.^2));
    dfcs = NSub*(NCond-1)
    MScs = SScs/dfcs;
    Fcont = MScont/MScs
    p = 1-fcdf(Fcont,1,dfcs)
end
