Find the correlation between 14 subjects

Hello,

I want to find the correlation between the 14 subjects that each subjects has 4 channels (columns). I wrote this 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.

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|')