跨频率耦合图解 |
您所在的位置：网站首页 › cfc滤波器 › 跨频率耦合图解 |
|
跨频率耦合图解¶
引言¶
先来通过一个交互式例子感受一下跨频率耦合,建立起一些直观的感受. 注解 点击下面的Activate,再点击run来启动交互 from ipywidgets import interactive import matplotlib.pyplot as plt from matplotlib.widgets import TextBox import numpy as np plt.rcParams['figure.figsize'] = (10.0, 10.0) def set_center_axis(ax): ax.spines['right'].set_color('none') ax.spines['top'].set_color('none') ax.spines['bottom'].set_position(('data', 0)) ax.spines['left'].set_position(('data', 0)) seconds = 1 sample_rate =500 def cfc_data(slow_ph=0,slow_ap=1,fast_ph=0,fast_ap=1): t = np.linspace(0, seconds, seconds*sample_rate) slow = np.sin(2*np.pi*5*t+slow_ph)*slow_ap fast = np.cos(2*np.pi*37*t+fast_ph)*fast_ap fast_am = 0.5*slow + 0.5 x = slow + fast * fast_am + np.random.randn(*t.shape) * .05 return slow,fast,fast_am,x,t def cfc_draw(slow_ph,slow_ap,fast_ph,fast_ap): slow,fast,fast_am,x,t=cfc_data(slow_ph,slow_ap,fast_ph,fast_ap) fig,(ax1,ax2,ax3,ax4)=plt.subplots(4,1) ax1.plot(t[:sample_rate], slow[:sample_rate]) textstr = '\n'.join((r'$slow=%.2f sin(2\cdot5\pi t+%.2f\pi)$' % (slow_ap,slow_ph/np.pi ),)) ax1.text(0.05,0.5,textstr,fontsize = 25) ax2.plot(t[:sample_rate], fast[:sample_rate]) textstr = '\n'.join((r'$fast=%.2f cos(2\cdot37\pi t+%.2f\pi)$' % (fast_ap,fast_ph/np.pi ),)) ax2.text(0.05,0.5,textstr,fontsize = 25) ax3.plot(t[:sample_rate], fast_am[:sample_rate]) ax3.text(0.05,0.5,'$fast\_am = 0.5\cdot slow + 0.5$',fontsize = 25) ax4.plot(t[:sample_rate], x[:sample_rate]) ax4.text(0.05,0.5,'$x=slow+fast\cdot fast\_am + 0.05random$',fontsize = 25) for ax in (ax1,ax2,ax3,ax4): set_center_axis(ax) def F(slow_ph,slow_ap,fast_ph,fast_ap): cfc_draw(slow_ph,slow_ap,fast_ph,fast_ap) interactive_plot = interactive(F,slow_ph=(0,2*np.pi,0.1*np.pi),slow_ap=(0,2,0.1),fast_ph=(0, 2*np.pi, 0.05*np.pi),fast_ap=(0,2,0.1)) interactive_plot跨频率耦合一共有三种形式: 相-相耦合（phase-phase coupling，PPC)； 相-幅耦合（phase-amplitude coupling，PAC）； 幅-幅耦合（amplitude-amplitude coupling，AAC）. 上面的交互式例子展示的是低频相位调制高频幅度,属于相-幅耦合. 如何对跨频率耦合现象进行分析和量化?¶目前研究得比较多的是相-幅耦合,其耦合程度可以用同步化指标\(PAC\)来计算: \[ P A C=\left|n^{-1} \sum_{t=1}^{n} a_{t} e^{i \phi_{t}}\right| \]其中\(n\)表示时间点的总数,\(t\)表示时间点,\(a_t\)表示高频信号在\(t\)时刻的功率,\(\phi_t\)表示低频信号在\(t\)时刻的相位,\(i\)为复数单位. 那么,跨频率耦合分析可以分三步走: 将信号分为高频段和低频段； 从滤波后的信号中提取振幅和相位； 确定相位和振幅是否相关(计算PAC). 步骤一:将信号分为高频段和低频段¶以上文交互式模组中的数据为例,进行分析. from pylab import * slow,fast,fast_am,x,t=cfc_data() figure(figsize=(14, 4)) # Create a figure with a specific size. plot(t[:sample_rate*1], x[:sample_rate*1]) []观察上图,可以猜测这个信号可能由一个5Hz的低频信号和一个25Hz以上的高频信号叠加而成. 对数据进行傅利叶变换,得到其功率谱,确认一下信号的频率成分. LFP = x dt = t[1] - t[0] # 定义采样间隔 T = t[-1] # ... 数据的持续时间, N = len(LFP) # ... 数据点数 x_h = hanning(N) * LFP # 将数据乘以汉宁窗 xf = rfft(x_h - x_h.mean()) # 计算傅里叶变换 Sxx = 2*dt**2/T * (xf*conj(xf)) # 计算功率谱 Sxx = real(Sxx) # 留下实部 df = 1 / T # 定义的频率分辨率, fNQ = 1 / dt / 2 # ... 和奈奎斯特频率 faxis = arange(0, fNQ + df, df) # 构造频率坐标 plt.plot(faxis, 10 * log10(Sxx))# 绘制频谱与频率. xlim([0, 50]) # 设置频率范围, ylim([-80, 0]) # ... 和功率范围 xlabel('Frequency [Hz]') # 给坐标轴打上标签 ylabel('Power [mV$^2$/Hz]');功率谱密度图显示,该信号主要由一个5Hz左右的低频和一个30~45Hz左右的高频组成.因此,信号可以划分为: 低频带: \([2,7]Hz\) 高频带: \([30,45]Hz\) from scipy import signal Wn = [2,7]; # 设置通频带[2-7]Hz， n = 100; # ... 和滤波器阶数, # ... 构建带通滤波器, b = signal.firwin(n, Wn, nyq=fNQ, pass_zero=False, window='hamming'); Vlo = signal.filtfilt(b, 1, LFP); # ... 将其应用到数据中. Wn = [30, 45]; # 设置通频带[30-45]Hz, n = 100; # ... 和滤波器阶数, # ... 构建带通滤波器, b = signal.firwin(n, Wn, nyq=fNQ, pass_zero=False, window='hamming'); Vhi = signal.filtfilt(b, 1, LFP); # ... 将其应用到数据中. figure(figsize=(14, 4)) plot(t, LFP) plot(t, Vlo) plot(t, Vhi) xlabel('Time [s]') xlim([24, 26]); ylim([-2, 2]); legend(['origin signal', 'low pass band signal', 'high pass band signal']); 步骤二:从滤波后的信号中提取振幅和相位¶在上一个步骤中,低频成分和高频成分被从原始信号中分解出来了,接下来提取低频成分的相位和高频成分的幅值. 在这个步骤中使用到了希尔伯特变换.关于希尔伯特变换的内容,笔者将在另一篇博客中展开. phi = angle(signal.hilbert(Vlo)) # 计算低频成分相位 amp = abs(signal.hilbert(Vhi)) # 计算高频成分的幅值 plot(amp) plot(Vhi) plot(phi) legend(['high freq signal envelope', 'high freq signal amplitude', 'low freq signal phase']); xlim([0,500]) (0.0, 500.0) 步骤三: 确定相位和振幅是否相关(计算PAC)¶ \[ P A C=\left|n^{-1} \sum_{t=1}^{n} a_{t} e^{i \phi_{t}}\right| \] p_bins = arange(-pi, pi, 0.1) a_mean = zeros(size(p_bins)-1) p_mean = zeros(size(p_bins)-1) for k in range(size(p_bins)-1): #对于每个相位区间, pL = p_bins[k] #... 左范围, pR = p_bins[k+1] #... 右范围. indices=(phi>=pL) & (phi=pL) & (phi h for s in hS]) / len(hS) print(p) 0.0 |
今日新闻 |
点击排行 |
|
推荐新闻 |
图片新闻 |
|
专题文章 |
| CopyRight 2018-2019 实验室设备网 版权所有 win10的实时保护怎么永久关闭 |