基于MI的cfc（交叉频率耦合）分析

PAC 指高节律振荡幅值 ( 能量 ) 受到低节律振荡的相位调制，被锁定在低节律振荡的相位上。PAC 现象是目前发现比较多、研究较为深入的一种CFC 现象。PAC 现象已在啮齿类动物的海马、基底核和猕猴的新皮质以及人类的皮层和海马中有发现。

MI 算法由Canolty Tort 等［70］发展提出，主要是基于高频节律的幅值在低频节律相位信号上的分布，来计算该分布的香农熵从而表征该分布的不均匀性，并且使用均匀分布时的香农熵最大值进行归一化．MI 值越大，代表低频节律的相位对高频节律的幅值的调节作用越强。Tort 等[72]同时采集并分析了大鼠在执行 T 迷宫任务时脑纹状体与海马处的 LFP信号，分析显示每处 LFP 均存在 theta-gamma 节律间 PAC ；并观察到大鼠执行迷宫任务过程中 PAC可迅速地从无耦合提升到很强耦合，再进一步恢复到基本无耦合状态，而在决策时 PAC 耦合强度达最大值。

具体的MI算法为：

通过希尔伯特变换提取时间序列x(t)的相位序列Φt

并用同样的方法提取时间序列y(t)的幅值包络序列At

将相位序列划分为n个等长的相位段，本文划分为18个相位段（从0°至360°），每段20，分别计算每个相位段上的序列的平均幅值，并将其中第j个相位段上分布的。

运用标准化熵度量H来表征调制指数，定义为

N为相位段的个数，pj定义为：

当时，可以得到H的最大可能值，通过标准化，可以得到基于熵的调制指数

相关函数

function [ mean_mi, max_mi , fre, MI] = PACMI( x,f1_pha,f2_pha,f1_mag,f2_mag ,step_w,windowlength)
%PACMI 求出信号的MI矩阵的最大值、平均值、以及最大值的相频和幅频的信息

% input---f1_pha=f1_pha:Start frequency of phase
% ---f2_pha=f2_pha;end frequency of phase
% ---f1_mag=f1_mag; Start frequency of magnitude
% ---f2_mag=f2_mag; end frequency of magnitude
% ---x: signal 16*m
% ----,step_w: step of window
% -----windowlength:length of window

% output ---mean_mi : the mean of MI
% --- max_mi : the max of MI
% --- fre : 2*16 ;1st row is magnitude frequency of max of MI
% --- 2ed row is phase frequency of max of MI
% ---MI : the MI
% 此处显示详细说明

f1_pha=f1_pha;
f2_pha=f2_pha;
f1_mag=f1_mag;
f2_mag=f2_mag;
fs=1000;

[ dataout ] = cfcfilt(x,fs,f1_pha,f2_mag,[],0 ,0 );

A=dataout;
[m,n]=size(A);

t = 1/fs:1/fs:n/fs;
node=f2_mag./step_w;
pac=zeros(2,n,node);

for i=1:node

LFP_all{i}=cfcfilt(A,fs,(i-1)*step_w, (i-1)*step_w+windowlength,[],0 ,0 );
for j=1:m
B=LFP_all{i}(j,:);
x = hilbert(B);
mag_x=abs(x);%求幅值，abs：求信号的模
pha_x=angle(x);%求相角
mag(j,:,i)=mag_x;
pha(j,:,i)=pha_x;
end
end
for ii=1:m % every channel
for jj=f1_pha:f2_pha %every phase band
for iii=f1_mag:f2_mag %every mag band
pac=[pha(ii,:,jj);mag(ii,:,iii)];
[y,H_p]=CFCMI(pac);
MI{ii}(iii,jj)=y;

HP{ii}{iii,jj}=H_p;
end
end
end
for jjj=1:m
pacmi(:,:)=MI{jjj};
mean_mi(1,jjj)=mean(mean(pacmi));
max_mi(1,jjj)=max(max(pacmi));
[x_f y_f]=find(pacmi==max(max(pacmi)));
fre(1,jjj)=x_f;
fre(2,jjj)=y_f;
end
end

function [ MI,H_p ] =CFCMI( x)
% CFCMI 将经过相幅信号x（2*m）（第一行为相位；第二行为幅值）
% 相位分为18段，每段20度，并求出个区间内所对应的幅值序列的平均值，
% 通过标准化熵度量H来表征调制指数，将Ｈ标准化得到MI
% input --x:2*m 相幅信号第一行为相位；第二行为幅值
% output --MI H_p 相幅耦合值
[m,n]=size(x);
pha_con=18;%18 phase band
step_pha=pi/pha_con;%Every band is 20 degrees
for i=1:pha_con
nar=find((x(1,:)>(i-1)*step_pha)&(x(1,:)<=i*step_pha));
pha_step{i}=x(:,nar);
A=pha_step{i};
B=A(2,:);
mag_mean=mean(B);
mag(:,i)=mag_mean;
end
p(1,:)=mag(1,:)/sum(mag);
for k=1:pha_con
cc=p(1,k);
H_p(1,k)=-cc.*log(cc);
end
H=sum(H_p);
pj=1/pha_con;
Hmax=-log(pj);%when pj=1/N, H is the max
MI=(Hmax-H)./Hmax;
end

function [ dataout ] = cfcfilt( data,fs,locutoff,hicutoff,filtorder, revfilt ,plotfreqz )
%UNTITLED5 此处显示有关此函数的摘要
%Inputs:
% data - data
% fs - frequency sample
% locutoff - lower edge of the frequency pass band (Hz)
% {[]/0 -> lowpass}
% hicutoff - higher edge of the frequency pass band (Hz)
% {[]/0 -> highpass}
%
% Optional inputs:
% filtorder - filter order (filter length - 1). Mandatory even
% revfilt - [0|1] invert filter (from bandpass to notch filter)
% {default 0 (bandpass)}
% usefft - ignored (backward compatibility only)
% plotfreqz - [0|1] plot filter's frequency and phase response
% {default 0}
% minphase - scalar boolean minimum-phase converted causal filter
% {default false}
%
% Outputs:
% dataout - filtered data
% 此处显示详细说明

[m,n]=size(data);

TRANSWIDTHRATIO = 0.25;
fNyquist = fs / 2;
edgeArray = sort([locutoff hicutoff]);
% Max stop-band width
maxTBWArray = edgeArray; % Band-/highpass
if revfilt == 0 % Band-/lowpass
maxTBWArray(end) = fNyquist - edgeArray(end);
elseif length(edgeArray) == 2 % Bandstop
maxTBWArray = diff(edgeArray) / 2;
end
maxTBWArray(end) = fNyquist - edgeArray(end);
maxDf = min(maxTBWArray);
if isempty(filtorder)
df = min([max([edgeArray(1) * TRANSWIDTHRATIO 2]) maxDf]);
filtorder = 3.3 / (df / fs); % Hamming window
filtorder = ceil(filtorder / 2) * 2; % Filter order must be even.
if filtorder==inf
filtorder=3.3*fs;
end

else

df = 3.3 / filtorder * fs; % Hamming window
filtorderMin = ceil(3.3 ./ ((maxDf * 2) / fs) / 2) * 2;
filtorderOpt = ceil(3.3 ./ (maxDf / fs) / 2) * 2;
end
filterTypeArray = {'lowpass', 'bandpass'; 'highpass', 'bandstop (notch)'};
fprintf('cfcfilt() - performing %d point %s filtering.\n', filtorder + 1, filterTypeArray{revfilt + 1, length(edgeArray)})
fprintf('cfcfilt() - transition band width: %.4g Hz\n', df)
fprintf('cfcfilt() - passband edge(s): %s Hz\n', mat2str(edgeArray))

% Passband edge to cutoff (transition band center; -6 dB)
dfArray = {df, [-df, df]; -df, [df, -df]};
cutoffArray = edgeArray + dfArray{revfilt + 1, length(edgeArray)} / 2;
fprintf('pop_eegfiltnew() - cutoff frequency(ies) (-6 dB): %s Hz\n', mat2str(cutoffArray))

% Window
winArray = windows('hamming', filtorder + 1);
b_f= cutoffArray / fNyquist;
% Filter coefficients
if b_f(1,1)==0
b_f(1,1)=1.000e-04;
end
b = firws(filtorder, b_f, winArray);

% Plot frequency response
if plotfreqz
freqz(b, 1, 8192, fs);
end

% Filter
disp('pop_eegfiltnew() - filtering the data (zero-phase)');

nFrames = 1000;
groupDelay = (length(b) - 1) / 2;
dcArray=[1,length(data)+1];

% Initialize progress indicator
nSteps = 20;
step = 0;
fprintf(1, 'firfilt(): |');
strLength = fprintf(1, [repmat(' ', 1, nSteps - step) '| 0%%']);
tic
for iDc = 1:(length(dcArray) - 1)
% Pad beginning of data with DC constant and get initial conditions
ziDataDur = min(groupDelay, dcArray(iDc + 1) - dcArray(iDc));
ans=double([data(:, ones(1, groupDelay) * dcArray(iDc)) ...
data(:, dcArray(iDc):(dcArray(iDc) + ziDataDur - 1))]);
[temp, zi] = filter(b, 1, ans, [], 2);

blockArray = [(dcArray(iDc) + groupDelay):nFrames:(dcArray(iDc + 1) - 1) dcArray(iDc + 1)];

for iBlock = 1:(length(blockArray) - 1)

% Filter the data
[data(:, (blockArray(iBlock) - groupDelay):(blockArray(iBlock + 1) - groupDelay - 1)), zi] = ...
filter(b, 1, double(data(:, blockArray(iBlock):(blockArray(iBlock + 1) - 1))), zi, 2);

% Update progress indicator
[step, strLength] = mywaitbar((blockArray(iBlock + 1) - groupDelay - 1), size(data, 2), step, nSteps, strLength);
end
temp = filter(b, 1, double(data(:, ones(1, groupDelay) * (dcArray(iDc + 1) - 1))), zi, 2);
data(:, (dcArray(iDc + 1) - ziDataDur):(dcArray(iDc + 1) - 1)) = ...
temp(:, (end - ziDataDur + 1):end);

% Update progress indicator
[step, strLength] = mywaitbar((dcArray(iDc + 1) - 1), size(data, 2), step, nSteps, strLength);

end
dataout=data;
end

基于MI的cfc（交叉频率耦合）分析相关推荐

workbench 流固耦合_基于Workbench的流固耦合作用下三通管振动特性分析
基于Workbench的流固耦合作用下三通管振动特性分析韩天宇,郭长青*,谌冉曦 (南华大学土木工程学院,湖南衡阳 421001) 摘要:使用ANSYS Workbench软件,对流固耦合作用 ...
abaqus dat文件 matlab_基于MPCCI的FLUENT与ABAQUS流固耦合分析步骤
FSI实例 FLUENT+ABAQUSMPCCI 以一个实例为例,说明如何采用多场耦合平台MPCCI将ANSYS Fluent的流场数据在每个计算步长内传递至ABAQUS中进行固体分析. 1.准备阶段 ...
解开神经科学中的交叉频率耦合
跨频率耦合(Cross-frequency coupling, CFC)已被提出在空间和时间尺度上协调神经动力学.尽管它对理解健康和病理的大脑功能有潜在的意义,但标准的CFC分析和生理学解释仍然存在根 ...
AI+医疗：基于模型的医疗应用大规模分析 | 腾讯AI Lab学术论坛演讲
来源:腾讯AI实验室摘要:3月15日,腾讯AI Lab第二届学术论坛在深圳举行,聚焦人工智能在医疗.游戏.多媒体内容.人机交互等四大领域的跨界研究与应用. 3月15日,腾讯AI Lab第二届学术论坛 ...
java消费者模式_基于Java 生产者消费者模式(详细分析)
生产者消费者模式是多线程中最为常见的模式:生产者线程(一个或多个)生成面包放进篮子里(集合或数组),同时,消费者线程(一个或多个)从篮子里(集合或数组)取出面包消耗.虽然它们任务不同,但处理的资源是相 ...
ansys流固耦合分析与工程实例_ansys workbench 流固耦合教程
点击蓝字关注我们流固耦合概念介绍流固耦合问题是流体力学(Computational Fluid Dynamics,CFD)与固体力学 (Computational Solid Mechanic ...
基于Android系统的IPv6网络接入分析
基于Android系统的IPv6网络接入分析摘要:本文深入分析了 ...
计算机模拟爆破过程,基于LSDYNA岩石爆破模拟建模分析
原标题:基于LSDYNA岩石爆破模拟建模分析 Ls Dyna岩石爆破模拟仿真作者:dyna_focus 擅长领域:dyna/abaqus/hypermesh 一. 数值模型的建立 1.1 单元及算法 ...
基于MATLAB的平面刚架有限元分析,基于MATLAB的平面刚架静力分析
工程计算实践作业基于MATLAB的平面刚架静力分析为了进一步理解有限元方法计算的过程,本文根据矩阵位移法的基本原理应用MATLAB编制计算程序对以平面刚架结构进行了静力分析.本文还利用ANSYS大 ...
基于BOOTSTRAP方法的投资方案分析
[摘要] 日常生活中存在各种各样的理财产品,本文以"华夏财富宝A"为例对常见经济现象进行数学建模并加以分析.实际情况下,基金收益情况受各种因素影响,但总体上遵循一定数学规律.本建 ...

基于MI的cfc（交叉频率耦合）分析

基于MI的cfc（交叉频率耦合）分析相关推荐

最新文章

热门文章