spike lfp.zip

function LFP_spike = lfp_spike_syn(bandpass_lfp,spk_times_stim,s_rate,dp_corr) % 输入： % bandpass_lfp：感兴趣频段的LFP信号； % spk_times_stim：spike的发放时刻； % s_rate：采样率； % dp_corr：数据点数，窗口长度*s_rate/2。 % 输出： % LFP_spike：LFP和spike的耦合强度。 N = 100; % 数据随机化次数 phase_ori =unwrap(angle(hilbert(bandpass_lfp))); M = length(spk_times_stim); for kk = 1:M spk_tmp = spk_times_stim(kk); phase_mat(kk,:) = phase_ori(round(spk_tmp*s_rate)-dp_corr:round(spk_tmp*s_rate)+dp_corr); %%% 初始相位矩阵 end DP = size(phase_mat,2); [C eig eig_vector] = C_phase_syn(phase_mat); for yy = 1:N for kk = 1:M idx1 = randperm(DP); tmp1 = phase_mat(kk,:); phase_surr(yy,kk,:) = tmp1(idx1); end end for yy = 1:N tmp = reshape(phase_surr(yy,:,:),M,DP); [C_surr eig_surr(yy,:)] = C_phase_syn(tmp); end mean_eig_surr = mean(eig_surr);%N次随机化后特征值的平均值 SD_eig_surr = std(eig_surr);%N次随机化后特征值的标准差 EG = eig - mean_eig_surr -3*SD_eig_surr;%特征值减去随机化的特征值 sumC=0; sumR=0; for s = M:-1:1 if EG(s)>0 sumC=sumC+EG(s);%所有大的特征值减去随机特征值的和 sumR=sumR+mean_eig_surr(s);%大的随机特征值的和 else break end end global_C = (eig - mean_eig_surr)./(M - mean_eig_surr); LFP_spike = global_C(end);