Find the correlation between 14 subjects

7 views (last 30 days)
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|')

Answers (0)

Categories

Find more on Vibration Analysis in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!