% demonstrates the use of phase slope index (PSI) as formulated in the paper:
%    Nolte G, Ziehe A, Nikulin VV, Schl\"ogl A, Kr\"amer N, Brismar T,
%    M\"uller KR. 
%    Robustly estimating the flow direction of information in complex
%    physical systems. 
%    Physical Review Letters. To appear. 
%    (for further information see    http://doc.ml.tu-berlin.de/causality/
%    )
%

% License
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see http://www.gnu.org/licenses/.







% trivial example for flow from channel 1 to channel 2. 
n=10000;x=randn(n+1,1);data=[x(2:n+1),x(1:n)]; 

% parameters for PSI-calculation 
segleng=100;epleng=200;

% calculation of PSI. The last argument is empty - meaning that 
% PSI is calculated over all frequencies
[psi, stdpsi, psisum, stdpsisum]=data2psi(data,segleng,epleng,[]);

% note, psi, as calculated by data2psi corresponds to \hat{\PSI} 
% in the paper, i.e., it is not normalized. The final is 
% the normalized version given by: 
psi./(stdpsi+eps)


% For psisum and stdpsisum refer to the paper


% To calculate psi in a band set, e.g., 
freqs=[5:10];
[psi, stdpsi, psisum, stdpsisum]=data2psi(data,segleng,epleng,freqs);
%with result:
psi./(stdpsi+eps)
% In this example, the flow is estimated to go from channel 
% 1 to channel 2 because the matrix element psi(1,2) is positive. 

% You can also calculate many bands at once, e.g. 
freqs=[[5:10];[6:11];[7:12]];
[psi, stdpsi, psisum, stdpsisum]=data2psi(data,segleng,epleng,freqs);
%and psi has then 3 indices with the last one refering to the row in freqs: 
psi./(stdpsi+eps)