%% THEORETICAL "BLOCK" EXAMPLE 
% references:
% L.Faes and G. Nollo, "Measuring Frequency Domain Granger Causality for Multiple Blocks of Interacting Time Series", Biological Cybernetics 2013. DOI: 10.1007/s00422-013-0547-5
% (this code realizes Fig.2 of the paper)
% L.Faes, S. Erla and G. Nollo, "Block Partial Directed Coherence: a New Tool for the Structural Analysis of Brain Networks", International Journal of Bioelectromagnetism, Vol. 14, No. 4, pp. 162 - 166, 2012


clear; close all; clc;
nfft=512; % number of frequency bins
fc=1; % sample frequency
simunumber=1;
Q=8; % number of series
Mv=[2 2 3 1]'; % number of series in each block
% % Mv=[1 1 1 1 1 1 1 1]'; % number of series in each block
pmax=2; % maximum lag

%% MVAR coefficients
ca=[0.5 0.5];
cb=[0.4 0.4];
cc=[0.4 0.4 0.4];
cd=[0.7];
ck=[-0.5 0.5 0.5 0.5];

pmax=2; % maximum lag
r1=0.9; f1=0.1; % oscillation channel 1
r3=0.9; f3=0.25; % oscillation channel 3

a1=ca(1); a2=ca(2);
b1=cb(1); b2=cb(2);
c1=cc(1); c2=cc(2); c3=cc(3);
d1=cd(1);
k1=ck(1); k2=ck(2); k3=ck(3); k4=ck(4);
Su=eye(Q); %Residual covariance matrix (DIAGONAL)
B0=zeros(Q,Q); %Matrix of instantaneous effects:no instantaneous effects
Bk=NaN*ones(Q,Q,pmax);  
%effects at lag 1
Bk(1,:,1)=[2*r1*cos(2*pi*f1) 0 0 0 0 0 0 0];
Bk(2,:,1)=[k1 0 0 0 0 0 0 0];
Bk(3,:,1)=[0 0 2*r3*cos(2*pi*f3) 0 0 0 0 d1];
Bk(4,:,1)=[0 a2 0 0 0 0 0 0];
Bk(5,:,1)=[0 0 c1 0 0 0 0 0];
Bk(6,:,1)=[0 0 0 0 k2 0 0 0];
Bk(7,:,1)=[0 0 c3 0 0 k4 0 0];
Bk(8,:,1)=[0 0 0 b2 0 0 0 0];
%effects at lag 2
Bk(1,:,2)=[-r1^2 0 0 0 0 0 0 0];
Bk(2,:,2)=[0 0 0 0 0 0 0 0];
Bk(3,:,2)=[0 a1 -r3^2 0 0 0 0 0];
Bk(4,:,2)=[0 0 0 0 0 0 0 0];
Bk(5,:,2)=[0 0 0 0 0 k3 0 0];
Bk(6,:,2)=[0 0 c2 0 0 0 0 0];
Bk(7,:,2)=[0 0 0 0 0 0 0 0];
Bk(8,:,2)=[0 0 b1 0 0 0 0 0];
% concateno in matrice Bm
Bm=[];
for kk=1:pmax
    Bm=[Bm Bk(:,:,kk)];
end
Am=Bm; % no inst eff
E=eye(Q*pmax);AA=[Am;E(1:end-Q,:)];lambda=eig(AA);
lambdamax=max(abs(lambda));

%% tmp: show one realization of the process
M=length(Mv);
Q=size(Am,1);
p=size(Am,2)/Q;

N=300;
U=randn(Q,N); % uncorrelated gaussian innovations
[Y]=MVARfilter(Am,U); % realization of Eq. 18
figure(1); % realizations of functions
for i=1:Q
    subplot(Q,1,i); plot(Y(i,:));
    xlim([0 N]);
    xlabel(['n']);ylabel(['x_' int2str(i)]);
end


%% Theoretical spectral functions
%%% block functions
[bDC,bPDC,mF,mG,bS,bP,bH,bAf,f] = block_fdMVAR(Am,Su,Mv,nfft,fc);
%[bDC2,bPDC2,mF2,mG2,bS,bP,bH,bAf,f] = block_fdMVAR_diag(Am,Su,Mv,nfft,fc);


%% graphs
h1=figure('numbertitle','off','name','Fig. 2a,b - Multivariate total and direct causality: f,g');
h2=figure('numbertitle','off','name','Fig. 2c,d - Block Directed Coherence, Block Partial Directed Coherence');

k=0;
for i=1:M        
    for j=i+1:M             
        figure(h1); % mF
        subplot(6,4,4*k+1); plot(f, squeeze(mF(i,j,:)),'b'); %hold on; plot(f, squeeze(mF2(i,j,:)),'c:');
        axis([0 fc/2 0 1.1*max(mF(:))]); title(['f^(^m^)' int2str(j) '->' int2str(i)]);
        subplot(6,4,4*k+2); plot(f, squeeze(mF(j,i,:)),'b'); %hold on; plot(f, squeeze(mF2(j,i,:)),'c:');
        axis([0 fc/2 0 1.1*max(mF(:))]); title(['f^(^m^)' int2str(i) '->' int2str(j)]);

        subplot(6,4,4*k+3); plot(f, squeeze(mG(i,j,:)),'r'); %hold on; plot(f, squeeze(mG2(i,j,:)),'m:');
        axis([0 fc/2 0 1.1*max(mF(:))]); title(['g^(^m^)' int2str(j) '->' int2str(i)]);
        subplot(6,4,4*k+4); plot(f, squeeze(mG(j,i,:)),'r'); %hold on; plot(f, squeeze(mG2(j,i,:)),'m:');
        axis([0 fc/2 0 1.1*max(mF(:))]); title(['g^(^m^)' int2str(i) '->' int2str(j)]);
        
        figure(h2); % mF
        subplot(6,4,4*k+1); plot(f, squeeze(bDC(i,j,:)),'b'); %hold on; plot(f, squeeze(bDC2(i,j,:)),'c:');
        axis([0 fc/2 -0.05 1.05]); title(['bDC' int2str(i) int2str(j)]);
        subplot(6,4,4*k+2); plot(f, squeeze(bDC(j,i,:)),'b'); %hold on; plot(f, squeeze(bDC2(j,i,:)),'c:');
        axis([0 fc/2 -0.05 1.05]); title(['bDC' int2str(j) int2str(i)]);

        subplot(6,4,4*k+3); plot(f, squeeze(bPDC(i,j,:)),'r'); %hold on; plot(f, squeeze(bPDC2(i,j,:)),'m:');
        axis([0 fc/2 -0.05 1.05]); title(['bPDC' int2str(i) int2str(j)]);
        subplot(6,4,4*k+4); plot(f, squeeze(bPDC(j,i,:)),'r'); %hold on; plot(f, squeeze(bPDC2(j,i,:)),'m:');
        axis([0 fc/2 -0.05 1.05]); title(['bPDC' int2str(j) int2str(i)]);

        k=k+1;        

    end
end
