MATLAB demo script for granger causality (powered by MVGC toolbox)

gc_spectral_demo.m

% 
% Jeelab Spectral Granger-causality demo script for dummies
% Toolbox source: http://www.sussex.ac.uk/sackler/mvgc/
% 
% written by Hio-Been Han, 2021-07-23. hiobeen.han@kaist.ac.kr
% https://jeelab.net
% 

%% (1) Environment setting
clear all;
addpath('mvgc_v1.0/'); startup;

%% (2) Demonstration parameters
actual_lag = 100; % unit: data point (regardless of srate)
signal_freq = [40, 60]; % informative frequency band (range, Hz)
noise_intensity = 1; % Signal intensity = approximately 1
lags_to_test = round([actual_lag*4, actual_lag*2, actual_lag, actual_lag/2, actual_lag/4]);

%% (3) Signal synthesis
srate = 1024;
signal_length = 3; % sec
N = srate*signal_length;
t = 1/srate:1/srate:N/srate;

rng(2021-07-23);
original_source = zerofilt( randn( [1, N+actual_lag] ), ...
    signal_freq(1), signal_freq(2), srate ) * 10;
x_source = original_source( actual_lag+1:N+actual_lag );
x_leader = original_source(actual_lag+1:N+actual_lag)+ randn([1,N]) * noise_intensity;
y_follower = original_source(1:N) + randn([1,N]) * noise_intensity;

%% (4) Calculation & visualization
tic
fig=figure(0723); clf; set(fig, 'Color', [1 1 1]); 
fig.Position(1)=300;fig.Position(2)=100;
fig.Position(3)=fig.Position(1)+1080; fig.Position(4)=fig.Position(2)+640;

% Signal visualization
subplot(2,1,1);
plot_multichan( t*1000, [x_source; x_leader; y_follower], 15 );
title( sprintf( 'Synthesized signal (signal=[%d-%d] Hz, actual lag=%03d (%02d ms in %d Hz), snr=%.1f)',...
    signal_freq(1), signal_freq(2), actual_lag, round(1000*actual_lag/srate), srate, 1/noise_intensity))
xlabel('Time (ms)');
set(gca, 'FontSize', 13,'Box', 'off', 'LineWidth', 2 );
drawnow;

% GCCA toolbox formatting
U = [ x_leader; y_follower ];

% Main calculation loop
for lagIdx = 1:length(lags_to_test)
    
    % Calculation
    lag_to_consider = lags_to_test(lagIdx);
    gc_xtoy = GCCA_tsdata_to_smvgc(U,2,1,lag_to_consider,srate/2-1);
    gc_ytox = GCCA_tsdata_to_smvgc(U,1,2,lag_to_consider,srate/2-1);
    freq_axis = linspace(0, srate/2, length(gc_xtoy));
    
    % Visualization
    subplot(2,length(lags_to_test),length(lags_to_test)+lagIdx); hold off;
    plt1=plot( freq_axis, gc_xtoy, 'LineWidth', 2 );
    hold on;
    plt2=plot( freq_axis, gc_ytox , 'LineWidth', 2);
    
    % Enhancing visibility
    xlabel('Frequency (Hz)');
    ylabel('Granger causality (a. u.)');
    xlim([1 100]);
    legend([plt1, plt2], 'X->Y (fwd)', 'Y->X (inv)', 'FontSize', 13);
    set(gca, 'FontSize', 13,'Box', 'off', 'LineWidth', 2 );
    title( sprintf( 'Considered lag = %03d', lag_to_consider ));
    
    % Setting uniform ylim to compare
    if lagIdx==1, ylims=ylim;  end
    ylim([ylims(1) ceil(ylims(2)*1.5)]);
    drawnow ;
end

toc

if true, saveas(gcf, 'gc_demo_figure.png'); end

%% Subfunctions
function result = zerofilt(data,hpfreq,lpfreq,srate)
Nyq=srate/2; f=[ 0 hpfreq hpfreq lpfreq lpfreq Nyq]/Nyq;  
a= [ 0 0 1 1 0 0];  b = firls(500,f,a); result = filtfilt(b,1,data);end
function [X,freq]=hiobeen_fft(x,srate)
N=(length(x));  k=0:N-1; T=N/srate; freq=k/T; cutOff = ceil(N/2); 
freq = freq(1:cutOff); X=fft(x)/N*2; X = X(1:cutOff);  end
function plot_multichan( t, x, interval )
nChan = size(x,1); stds = nanstd( x, 0, 2 );for chIdx = 1:size(x,1), x(chIdx,:) = nanmean(stds) * (x(chIdx,:) / stds(chIdx)); end
y_center = linspace( -interval, interval, nChan );
c_space=imresize(colormap('lines'),[nChan,3],'nearest');
chanlab = {['Y: Follow+' char(949)], ['X: Lead+' char(949)], 'Signal source' }; chanlab_pos = [];
for chanIdx = 1:nChan
    plot( t, x( chanIdx, : ) - (y_center(chanIdx) + nanmean( x( chanIdx, : ), 2)), 'Color', c_space( chanIdx,: )); 
    chanIdx_reverse = nChan-chanIdx+1; chanlab_pos(chanIdx) =  y_center(chanIdx) ; if chanIdx ==1, hold on; end
end; hold off; set(gca, 'YTick', chanlab_pos, 'YTickLabel', chanlab); ylim([-1 1]*interval*1.5); end