How to find the correlation between 14 subjects

Question

Neda Deljavan 2022-10-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1821480-how-to-find-the-correlation-between-14-subjects

评论： Benjamin Thompson 2022-10-17

Hello,

I want to find the correlation between the 14 subjects that each subjects has 4 channels (columns). I wrote its code. But I do not know it is correct or not. The main problem for me is that "how can I find the correlation between 14 subjects". I mean I want to find the first square for subject1 which is shown in the figure below, then upto subject12.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Benjamin Thompson 2022-10-10

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1821480-how-to-find-the-correlation-between-14-subjects#answer_1070590

xcorr2 is the MATLAB function for 2D cross correlation

6 个评论
显示 4更早的评论隐藏 4更早的评论

Neda Deljavan 2022-10-11

在 MATLAB Online 中打开

Actually I coded this one. but I do want to code for all groups.

Here, I found correlation for channel1 of the first subject then channel2 of the first subject then channel3 of the first subject then channel4 of the first subject in a single code for each of them like the code below:

%%
clc
clear
close all
%% Data Description
% 30800 timepoints, 256Hz(sampling rate), 120s
% 4 channels(Lateral Frontal: ch1=F3, ch2=F7, ch3=F4, ch4=F8)
% 14 Subject(Active) * 4 conditions(cd1=EC1, cd2=EC2, cd3=EO1, cd4=EO2)= 56 
%% Step1: Load data & Plot
% Taking one channel F3
% Taking first condition= EC1
fs=256; % Sampling frequency(Hz)
load P01EC1
sig= P01EC1(:,1)';
% sig1= 3*(sin(2*pi*(1:1:30800)/fs*5)) + 0*(sin(2*pi*(1:1:30800)/fs*12));
figure
% subplot(4,1,1)
plot(sig,'b','linewidth',0.5)
title('P01,EC1Active,Chennel F3 in Time domain')
xlabel('Time (s)')
ylabel('Voltage (mV)')
grid on
grid minor
%% Step2: Power Spectrum
% we plot a power to see peaks
N1 = length(sig);
fx1 = fft(sig);
fx1 = fx1(1:N1/2+1);
psd = (1/(fs*N1)) * abs(fx1).^2;
psd(2:end-1) = 2*psd(2:end-1);
freq = 0:fs/length(sig):fs/2;
figure
plot(freq,10*log10(psd))
xlim([0 55])
grid on
grid minor
title('P01,EC1Active,Periodogram Using FFT')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')
%% Step4: Denoising (50Hz)
% Step 4.1: Design stop filter
fl= 49.9;
fh= 50.1;
order= 3;
wn= [fl fh]/ (fs/2);
type= 'stop';
[b1,a1]= butter(order,wn,type);
%  Step 4.2: Apply notch (stop) filter
sig= filtfilt(b1,a1,sig);
%% Step5: Fourier Transform
N1= length(sig);
fx1= fft(sig,N1);
figure
% subplot(4,1,2)
stem(abs(fx1),'b','linewidth',0.5,'marker','none')
title('P01,EC1Active,Chennel F3 in Frequency domain')
grid on
grid minor
%% Step6: Select half of coeficients
fx1= fx1(1: round(N1/2));
figure
stem(abs(fx1),'b','linewidth',0.5,'marker','none')
title('P01,EC1Active,Chennel F3 in Frequency domain')
grid on
grid minor
%% Step 7 : Calculate Magnitude of coeficients
px1= abs(fx1);
%% Step 8 : Calculate Frequency Resolution
rf= linspace(0,fs/2,round(N1/2));
figure
stem(rf,px1,'b','linewidth',0.5,'marker','none')
xlim([0 55])
title('P01EC1,Active,Chennel F3 in Frequency domain')
grid on
grid minor
hold on
%% Step9: Determine coeficients of rhythms
% delta= 0.1-4 Hz,
% theta= 4-8 Hz,
% alpha= 8-12 Hz,
% beta= 12-20 Hz,
% lowgamma= 20-25 Hz,
% gamma= 25-40 Hz
name={'Delta','Theta','Alpha','Beta','Lowgamma','Gamma'};
band= [0.1,4,8,12,20,25;
    4,8,12,20,25,40];
figure
hold on
for j=1:size(band,2) %first column:delta
    fl= band(1,j);
    fh= band(2,j);
    indx1= find(rf>=fl & rf<fh);
    %         figure
    %         subplot(8,1,5)
    stem(rf(indx1),px1(indx1),'linewidth',0.5,'marker','none')
    title('Frequency Rhythms P01,EC1Active,Chennel F3 in Frequency domain')
    xpos=rf(indx1(1));
    ypos=max(px1(indx1));
    text(xpos,ypos,name{j},'fontsize',10)
    %         drawnow
end
hold off
%% Step10: Band decompotition
% Step 10.1: Plot the original signal
sig= P01EC1(:,1);
% figure
subplot(7,1,1)
plot(sig,'b','linewidth',0.5)
title('P01,EC1Active,Chennel F3 in Time domain')
xlabel('Time (s)')
ylabel('Voltage (mV)')
grid on
grid minor
% Step 10.2: Design Filter for delta
% bandpass butterworth filter
fl= 1;
fh= 4;
order= 3;
wn= [fl fh]/(fs/2);
type= 'bandpass';
[b,a]= butter(order,wn,type);
% Apply designed filter
sig1_delta= filtfilt(b,a,sig);
% Display result
% figure
subplot(7,1,2)
plot(sig1_delta,'b','linewidth',1)
title('Delta wave in Time domain')
grid on
grid minor
% Step 10.3: Design Filter for theta
% bandpass butterworth filter
fl= 4;
fh= 8;
order= 3;
wn= [fl fh]/(fs/2);
type= 'bandpass';
[b,a]= butter(order,wn,type);
% Apply designed filter
sig1_theta= filtfilt(b,a,sig);
% Display result
% figure
subplot(7,1,3)
plot(sig1_theta,'b','linewidth',1)
title('Theta wave in Time domain')
grid on
grid minor
% Step 10.4: Design Filter for alpha
% bandpass butterworth filter
fl= 8;
fh= 12;
order= 3;
wn= [fl fh]/(fs/2);
type= 'bandpass';
[b,a]= butter(order,wn,type);
% Apply designed filter
sig1_alpha= filtfilt(b,a,sig);
% Display result
% figure
subplot(7,1,4)
plot(sig1_alpha,'b','linewidth',1)
title('Alpha wave in Time domain')
grid on
grid minor
% Step 10.5: Design Filter for beta
% bandpass butterworth filter
fl= 12;
fh= 20;
order= 3;
wn= [fl fh]/(fs/2);
type= 'bandpass';
[b,a]= butter(order,wn,type);
% Apply designed filter
sig1_beta= filtfilt(b,a,sig);
% Display result
figure
subplot(7,1,5)
plot(sig1_beta,'b','linewidth',1)
title('Beta wave in Time domain')
grid on
grid minor
% Step 10.6: Design Filter for lowgamma
% bandpass butterworth filter
fl= 20;
fh= 25;
order= 3;
wn= [fl fh]/(fs/2);
type= 'bandpass';
[b,a]= butter(order,wn,type);
% Apply designed filter
sig1_lowgamma= filtfilt(b,a,sig);
% Display result
% figure
subplot(7,1,6)
plot(sig1_lowgamma,'b','linewidth',1)
title('Lowgamma wave in Time domain')
grid on
grid minor
% Step 10.7: Design Filter for gamma
% bandpass butterworth filter
fl= 25;
fh= 40;
order= 3;
wn= [fl fh]/(fs/2);
type= 'bandpass';
[b,a]= butter(order,wn,type);
% Apply designed filter
sig1_gamma= filtfilt(b,a,sig);
% Display result
% figure
subplot(7,1,7)
plot(sig1_gamma,'b','linewidth',1)
title('Gamma wave in Time domain')
grid on
grid minor
%% step11: Gathering all subjects in cd=1 EC1 Active
% % Load EC1, Active
% % B is the whole EC1 Active subjects
% 
% B=load('B.mat');
% sig= B(:,:);
% cd=1; % cd1= EC1 Active
% % subs_al= zeros(size(sig));
% % convert struct to double matrix
% B = struct2cell(sig);
% out = cat(2,B{:});
% % B= cell2mat(struct2cell(sig));
% Find Cov and Corr between one channel and all channels of one subject
load P01EC1
sig1= P01EC1(:,:);
%% step12: Design Filter for alpha
% % bandpass butterworth filter
% fl= 8;
% fh= 12;
% order= 3;
% wn= [fl fh]/(fs/2);
% type= 'bandpass';
% [b,a]= butter(order,wn,type);
% for i=1:14
%     ts_to_filt=sig.output_cell{i}; %cell to double matrix
%     subs_cd1_al=filtfilt(b,a,ts_to_filt);
% end
fl= 8;
fh= 12;
order= 3;
wn= [fl fh]/(fs/2);
type= 'bandpass';
[b,a]= butter(order,wn,type);
% for i=1:4
%     subs_EC1_alpha=filtfilt(b,a,sig1(:,i));
% end
subs_EC1_alpha=filtfilt(b,a,sig1);
%% step13: Covariance
cov_EC1 = cov(subs_EC1_alpha(:,:));
cov_EC1 = abs(cov_EC1 - diag(diag(cov_EC1)));
figure;
imagesc(cov_EC1);
axis square;
cb = colorbar;
cb.Label.String = 'Covariance level';
title('|Covariance, P01,EC1Active|')
ylabel('nodes');
xlabel('nodes')
%% step14: Pearson
% Linear correlaton between two random variables.
% In other words, how close the interdependency is near to a linear relation.
% It's also known to be the Covariance normalized as well as the $\tau=0$ of the Cross-Correlation function.
pea_EC1 = corrcoef(subs_EC1_alpha(:,:));
pea_EC1 = abs(pea_EC1 - diag(diag(pea_EC1)));
figure;
imagesc(pea_EC1);
axis square;
colorbar;
title('|Pearson, P01,EC1Active|')

Neda Deljavan 2022-10-12

编辑：Neda Deljavan 2022-10-12

The main question is that, How can I find the correlation between 14 subjects (each subject has 4channels as columns of matrix)?

In fact, I want "A matrix" that includes correlation between 14 subjects.

Benjamin Thompson 2022-10-17

Maybe you could post the data for two or three subjects and the Community could help you extend your code to do that part of the problem? For the A matrix construction, are you expecting the size to be 14*4 or 14? You can create an empty A matrix to fill with results using the zeros function:

A = zeros(14,14);

请先登录，再进行评论。