Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add more files, redundancy reduced #28

Merged
merged 18 commits into from
Mar 8, 2023
Merged

Add more files, redundancy reduced #28

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
18 commits
Select commit Hold shift + click to select a range
df3dc46
Update README.md
nguyen-td Jan 31, 2023
274aa14
Merge pull request #20 from arnodelorme/master
nguyen-td Feb 3, 2023
a992e91
add new test file, add docstring
nguyen-td Feb 3, 2023
0d68dd9
Merge branch 'feature' of https://github.com/Hiyeri/roiconnect into f…
nguyen-td Feb 3, 2023
5d1cc51
add test script for connectivity computation
nguyen-td Feb 3, 2023
a80c7a8
simplify code
nguyen-td Feb 3, 2023
5b80d39
Merge pull request #21 from arnodelorme/master
nguyen-td Feb 10, 2023
72390cf
Merge pull request #22 from Hiyeri/feature
nguyen-td Feb 10, 2023
9656fd4
fix wrong value to parameter list
nguyen-td Feb 22, 2023
a1b609d
Merge branch 'feature' of https://github.com/Hiyeri/roiconnect into f…
nguyen-td Feb 22, 2023
843ebe9
Merge pull request #23 from Hiyeri/dev
nguyen-td Feb 22, 2023
3aa28f1
- implemented wPLI in data2cs_event, needs to be tested first
nguyen-td Feb 22, 2023
902b3da
wPLI now computed in data2cs_event instead of in tsdata_to_cpsdwpli
nguyen-td Feb 24, 2023
fa2234f
Merge pull request #24 from Hiyeri/feature
nguyen-td Feb 24, 2023
18e826e
Merge pull request #25 from sccn/master
nguyen-td Feb 28, 2023
fd91a14
reduce redundant computation of cross-spectrum
nguyen-td Feb 28, 2023
819741d
Merge pull request #26 from Hiyeri/feature
nguyen-td Feb 28, 2023
ac63831
Merge pull request #27 from Hiyeri/dev
nguyen-td Mar 1, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
231 changes: 128 additions & 103 deletions libs/haufe/data2cs_event.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,36 +1,36 @@
function [cs,coh,nave]=data2cs_event(data,segleng,segshift,epleng,maxfreqbin,para);
% usage: [cs,coh,nave]=data2cs_event(data,segleng,segshift,epleng,maxfreqbin,para)
%
function [cs, coh, wpli, nave]=data2cs_event(data,segleng,segshift,epleng,maxfreqbin,para);
% usage: [cs, coh, wpli, nave]=data2cs_event(data,segleng,segshift,epleng,maxfreqbin,para)
%
% calculates cross-spectra from data for event-related measurement
% input:
% input:
% data: ndat times nchan matrix each colum is the time-series in one
% channel;
% segleng: length of each segment in bins, e.g. segleng=1000;
% segleng: length of each segment in bins, e.g. segleng=1000;
% segshift: numer of bins by which neighboring segments are shifted;
% e.g. segshift=segleng/2 makes overlapping segments
% epleng: length of each epoch
% maxfreqbin: max frequency in bins
% para: optional structure:
% para.segave=0 -> no averaging across segments
% para.segave=0 -> no averaging across segments
% para.segave neq 0 -> averaging across segments (default is 0)% \
% para.subave =1 subtracts the average across epochs,
% para.subave ~= 1 -> no subtraction (default is 1)
% para.subave =1 subtracts the average across epochs,
% para.subave ~= 1 -> no subtraction (default is 1)
% IMPORTANT: if you just one epoch (e.g. for continuous data)
% set para.subave=0
%
% set para.subave=0
%
% -> averaging across segments (default is 0)
% para.proj must be a set of vector in channel space,
% if it exists then the output raw contains the single trial
% Fourier-transform in that channel
%
%
% output:
% cs: nchan by chan by maxfreqbin by nseg tensor cs(:,:,f,i) contains
% para.proj must be a set of vector in channel space,
% if it exists then the output raw contains the single trial
% Fourier-transform in that channel
%
%
% output:
% cs: nchan by chan by maxfreqbin by nseg tensor cs(:,:,f,i) contains
% the cross-spectrum at frequency f and segment i
% coh: complex coherency calculated from cs
% coh: complex coherency calculated from cs
% nave: number of averages

subave=1;
subave=1;

if nargin<6
para=[];
Expand All @@ -41,47 +41,47 @@
segave=0;
mydetrend=0;
proj=[];
if isfield(para,'segave')
if isfield(para,'segave')
segave=para.segave;
end
if isfield(para,'detrend')
end
if isfield(para,'detrend')
mydetrend=para.detrend;
end
if isfield(para,'proj')
end
if isfield(para,'proj')
proj=para.proj;
end
if isfield(para,'subave')
end
if isfield(para,'subave')
subave=para.subave;
end
end

[ndum,npat]=size(proj);

[ndat,nchan]=size(data);
if npat>0
data=data*proj;
nchan=npat;
if npat>0
data=data*proj;
nchan=npat;
end

nep=floor(ndat/epleng);

nseg=floor((epleng-segleng)/segshift)+1; %total number of segments

if segave==0
cs=zeros(nchan,nchan,maxfreqbin,nseg);
av=zeros(nchan,maxfreqbin,nseg);
cs=zeros(nchan,nchan,maxfreqbin,nseg);
av=zeros(nchan,maxfreqbin,nseg);
else
cs=zeros(nchan,nchan,maxfreqbin);
av=zeros(nchan,maxfreqbin);
cs=zeros(nchan,nchan,maxfreqbin);
av=zeros(nchan,maxfreqbin);
end

if npat>0
if segave==0
cs=zeros(nchan,nchan,maxfreqbin,nep,nseg);
av=zeros(nchan,maxfreqbin,nep,nseg);
else
cs=zeros(nchan,nchan,maxfreqbin,nep);
av=zeros(nchan,maxfreqbin,nep);
end
if segave==0
cs=zeros(nchan,nchan,maxfreqbin,nep,nseg);
av=zeros(nchan,maxfreqbin,nep,nseg);
else
cs=zeros(nchan,nchan,maxfreqbin,nep);
av=zeros(nchan,maxfreqbin,nep);
end
end

mywindow=repmat(hanning(segleng),1,nchan);
Expand All @@ -90,96 +90,121 @@
end

%figure;plot(mywindow);
nave=0;
for j=1:nep;
nave = 0;
wpli_numer = cs; wpli_denom = cs; % initialize in the same way as CS

for j=1:nep
dataep=data((j-1)*epleng+1:j*epleng,:);
for i=1:nseg; %average over all segments;
for i=1:nseg %average over all segments;
dataloc=dataep((i-1)*segshift+1:(i-1)*segshift+segleng,:);
if mydetrend==1
datalocfft=fft(detrend(dataloc,0).*mywindow);
datalocfft=fft(detrend(dataloc,0).*mywindow);
else
datalocfft=fft(dataloc.*mywindow);
datalocfft=fft(dataloc.*mywindow);
end

for f=1:maxfreqbin % for all frequencies
if npat==0
if segave==0
cs(:,:,f,i)=cs(:,:,f,i)+conj(datalocfft(f,:)'*datalocfft(f,:));
av(:,f,i)=av(:,f,i)+conj(datalocfft(f,:)');
else
%disp([i,f,size(datalocfft)])
cs(:,:,f)=cs(:,:,f)+conj(datalocfft(f,:)'*datalocfft(f,:));
av(:,f)=av(:,f)+conj(datalocfft(f,:)');
end
else
if segave==0
cs(:,:,f,j,i)=conj(datalocfft(f,:)'*datalocfft(f,:));
av(:,f,j,i)=conj(datalocfft(f,:)');
else
%disp([i,f,size(datalocfft)])
cs(:,:,f,j)=cs(:,:,f,j)+conj(datalocfft(f,:)'*datalocfft(f,:));
av(:,f,j)=av(:,f,j)+conj(datalocfft(f,:)');
end
end

for f=1:maxfreqbin % for all frequencies
if npat==0
if segave==0
cs(:,:,f,i) = cs(:,:,f,i) + conj(datalocfft(f,:)'*datalocfft(f,:));
av(:,f,i) = av(:,f,i) + conj(datalocfft(f,:)');
wpli_numer(:,:,f,i) = wpli_numer(:,:,f,i) + imag(conj(datalocfft(f,:)'*datalocfft(f,:)));
wpli_denom(:,:,f,i) = wpli_denom(:,:,f,i) + abs(imag(conj(datalocfft(f,:)'*datalocfft(f,:))));
else
%disp([i,f,size(datalocfft)])
cs(:,:,f) = cs(:,:,f) + conj(datalocfft(f,:)'*datalocfft(f,:));
av(:,f) = av(:,f) + conj(datalocfft(f,:)');
wpli_numer(:,:,f) = wpli_numer(:,:,f) + imag(conj(datalocfft(f,:)'*datalocfft(f,:)));
wpli_denom(:,:,f) = wpli_denom(:,:,f) + abs(imag(conj(datalocfft(f,:)'*datalocfft(f,:))));
end
else
if segave==0
cs(:,:,f,j,i) = conj(datalocfft(f,:)'*datalocfft(f,:));
av(:,f,j,i) = conj(datalocfft(f,:)');
wpli_numer(:,:,f,j,i) = imag(conj(datalocfft(f,:)'*datalocfft(f,:)));
wpli_denom(:,:,f,j,i) = abs(imag(conj(datalocfft(f,:)'*datalocfft(f,:))));
else
%disp([i,f,size(datalocfft)])
cs(:,:,f,j) = cs(:,:,f,j) + conj(datalocfft(f,:)'*datalocfft(f,:));
av(:,f,j) = av(:,f,j)+conj(datalocfft(f,:)');
wpli_numer(:,:,f,j) = wpli_numer(:,:,f,j) + imag(conj(datalocfft(f,:)'*datalocfft(f,:)));
wpli_denom(:,:,f,j) = wpli_denom(:,:,f,j) + abs(imag(conj(datalocfft(f,:)'*datalocfft(f,:))));
end
end
end
end
nave=nave+1;
end

if segave==0
cs=cs/nave;
av=av/nave;
cs = cs/nave;
av = av/nave;
wpli_numer = wpli_numer/nave;
wpli_denom = wpli_denom/nave;
else
nave=nave*nseg;
cs=cs/nave;
av=av/nave;
nave = nave*nseg;
cs = cs/nave;
av = av/nave;
wpli_numer = wpli_numer/nave;
wpli_denom = wpli_denom/nave;
end

% compute wPLI
wpli = wpli_numer ./ wpli_denom;
wpli = permute(wpli, [3, 1, 2]);

for f=1:maxfreqbin
if subave==1
if npat==0
if segave==0
for i=1:nseg;cs(:,:,f,i)=cs(:,:,f,i)-av(:,f,i)*av(:,f,i)';end;
else
cs(:,:,f)=cs(:,:,f)-av(:,f)*av(:,f)';
end
else
if segave==0
for i=1:nseg;for j=1:nep;
cs(:,:,f,j,i)=cs(:,:,f,j,i)-av(:,f,j,i)*av(:,f,j,i)';
end;end;
else
for j=1:nep;cs(:,:,f,j)=cs(:,:,f,j)-av(:,f,j)*av(:,f,j)';end
end
end
end
if subave==1
if npat==0
if segave==0
for i=1:nseg
cs(:,:,f,i)=cs(:,:,f,i)-av(:,f,i)*av(:,f,i)';
end
else
cs(:,:,f)=cs(:,:,f)-av(:,f)*av(:,f)';
end
else
if segave==0
for i=1:nseg
for j=1:nep
cs(:,:,f,j,i)=cs(:,:,f,j,i)-av(:,f,j,i)*av(:,f,j,i)';
end
end
else
for j=1:nep
cs(:,:,f,j)=cs(:,:,f,j)-av(:,f,j)*av(:,f,j)';
end
end
end
end
end

ndim=length(size(cs));
if ndim==3;
if ndim==3
[n1,n2,n3]=size(cs);
coh=cs;
for i=1:n3;
for i=1:n3
c=squeeze(cs(:,:,i));
coh(:,:,i)=c./sqrt(diag(c)*diag(c)');
end
elseif ndim==4
[n1,n2,n3,n4]=size(cs);
coh=cs;
for i=1:n3;for j=1:n4;
c=squeeze(cs(:,:,i,j));
coh(:,:,i,j)=c./sqrt(diag(c)*diag(c)');
end;end;
for i=1:n3
for j=1:n4
c=squeeze(cs(:,:,i,j));
coh(:,:,i,j)=c./sqrt(diag(c)*diag(c)');
end
end
end








return;


7 changes: 3 additions & 4 deletions libs/haufe/data2sctrgcmim.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -100,15 +100,14 @@
CSpara.mywindow = hanning(ndat)./sqrt(hanning(ndat)'*hanning(ndat));


clear TRGC GC MIM MIC CS COH
clear TRGC GC MIM MIC CS COH wPLI

if abs(nboot) < 1 % no bootstrap

% data to autocovariance
% G = tsdata_to_autocov(data, nlags);
% CS = tsdata_to_cpsd_fast(data, fres, 'WELCH', ndat);
CS = data2cs_event(data(:, :)', ndat, floor(ndat/2), ndat, [], CSpara);

[CS, ~, wPLI, ~] = data2cs_event(data(:, :)', ndat, floor(ndat/2), ndat, [], CSpara);

maxfreq = size(CS,3);

if ~isempty(intersect(output, {'MIM', 'MIC', 'COH'}))
Expand Down
20 changes: 4 additions & 16 deletions libs/haufe/data2spwctrgc.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,7 @@
end

freqs = freqs(1:maxfreq);
z = exp(-i*pi*freqs);
z = exp(-1i*pi*freqs);

if nlags < 0
if cond
Expand Down Expand Up @@ -105,21 +105,9 @@

if verbose; disp(['computing cross-spectrum']); end

% crossspectrum using multitapers
if ~isempty(intersect(output, {'wPLI'}))
[CS wPLI] = tsdata_to_cpsdwpli(data, fres, 'MT', ndat, [], 1);
% [CS wPLI] = tsdata_to_cpsdwpli(data, fres, 'MT', (floor(ndat/2)), [], 1);

maxfreq = size(CS,3);
wPLI = wPLI(:, :, 1:maxfreq);
wPLI = permute(wPLI, [3, 1, 2]);
wPLI(isnan(wPLI)) = 0;
else
% CS = tsdata_to_cpsd(data, fres, 'MT', ndat, [], 1);
CS = data2cs_event(data(:, :)',ndat, ndat - (floor(ndat/2)), ndat , [], CSpara);
ASR = log(norm(imag(CS(:)))/norm(real(CS(:))));
fprintf('ASR: %1.4f\n', ASR);
end
CS = data2cs_event(data(:, :)',ndat, ndat - (floor(ndat/2)), ndat , [], CSpara);
ASR = log(norm(imag(CS(:)))/norm(real(CS(:))));
fprintf('ASR: %1.4f\n', ASR);

if ~isempty(intersect(output, {'COH'}))
clear COH
Expand Down
Loading