一种脑磁信号时频分布信息求解方法

1.本发明涉及分布信息求解领域，具体涉及一种脑磁信号时频分布信息求解方法。


背景技术：

2.脑磁信号属于多分量非稳态信号，信号频率活动范围广而且时变，同时大脑振荡的频率和出现的时间和大脑状态息息相关。随着硬件和预处理技术的发展，我们能够越来越准确地获得信噪比较高的脑磁信号，因此如何从信号中更准确地提取脑磁信号在时间-频率维度的特征，是脑磁系统走向应用的关键。其中，时频分析技术是将脑磁信号从时域映射到时-频域并计算信号激活情况的关键。目前，各个脑磁分析工具包及研究人员们常用的时频分析方法及其优缺点为：(1)stft：计算快，但是受固定的窗口大小的影响，时-频(tf)能量不集中，各个频段能量集中较差。(2)wt：小波在高频段具有良好的时间分辨率，但频率分辨率较差，在低频段则相反。该方法无法同时兼顾时间和频率分辨率。


技术实现要素：

3.针对现有技术的不足，本发明提供了一种脑磁信号时频分布信息求解方法，通过精细线性调频波估计，更加准确地对脑磁信号进行信号分离，准确提供了有关脑磁信号的时间频率的能量分布信息，能够同时较高分辨率地得到高频和低频时间分布特征。
4.为实现上述目的，本发明提供了一种脑磁信号时频分布信息求解方法，包括：
5.s1、获取初始脑磁信号进行降噪处理得到基础脑磁信号；
6.s2、利用所述基础脑磁信号进行分解处理得到脑磁啁啾波变换信号；
7.s3、利用所述脑磁啁啾波变换信号得到基于啁啾波变换的脑磁啁啾波变换数据时频表征；
8.s4、利用所述脑磁啁啾波变换数据时频表征得到脑磁信号瞬时激活频率。
9.优选的，所述获取初始脑磁信号进行降噪处理得到基础脑磁信号包括：
10.获取初始脑磁信号进行降噪处理得到的基础脑磁信号的计算式如下：
[0011][0012]
其中，x(t)为基础脑磁信号，t为时间，xk(t)是组成的第k个成分，n为成分个数，amk(t)为x(t)第k个成分的瞬时振幅，为x(t)第k个成分的瞬时相位，e为常数。
[0013]
优选的，利用所述基础脑磁信号进行分解处理得到脑磁啁啾波变换信号包括：
[0014]
s2-1、利用所述基础脑磁信号基于mp算法得到脑磁啁啾波变换粗估计信号；
[0015]
s2-2、利用所述脑磁啁啾波变换粗估计信号基于em算法得到脑磁啁啾波变换精估计信号；
[0016]
s2-3、利用所述脑磁啁啾波变换精估计信号得到脑磁啁啾波变换信号。
[0017]
进一步的，利用所述基础脑磁信号基于mp算法得到脑磁啁啾波变换粗估计信号包
括：
[0018]
利用所述基础脑磁信号计算基础脑磁信号的残差信号权重因子的计算式如下：
[0019][0020]
其中，ai为残差信号权重因子，xi(t)为第i个基础脑磁信号，gi(t)为第i个啁啾变换的啁啾波，g
i*
(t)为gi(t)的共轭复数；
[0021]
利用所述残差信号权重因子计算基础脑磁信号的残差信号的计算式如下：
[0022]ri+1
(t)＝xi(t)-a
igi
(t)
[0023]
其中，r
i+1
(t)为啁啾变换后的结果相对xi(t)的残差，xi(t)为第i个变换的基础脑磁信号成分，ai为残差信号权重因子，gi(t)为第i个啁啾变换的啁啾波；
[0024]
利用所述残差信号计算残差信号能量的计算式如下：
[0025][0026]
其中，r
i+1
(t)为残差信号能量，r
i+1*
(t)为r
i+1
(t)的共轭复数；
[0027]
其中，用于脑磁啁啾波变换的啁啾波基的计算式如下：
[0028][0029]
其中，g(t)为脑磁啁啾波变换的啁啾波基，tc为啁啾波时间中心，fc为啁啾波的频率中心，δt为采样间隔时间，c为啁啾率，d与δt相同，t为时间，j为虚数的虚部，i为虚数；
[0030]
利用所述脑磁啁啾波变换的参数基于mp算法得到脑磁啁啾波变换粗估计信号；
[0031]
其中，脑磁啁啾波变换的估计参数为tc、fc、log(δt)、c、d。
[0032]
进一步的，利用所述脑磁啁啾波变换粗估计信号基于em算法得到脑磁啁啾波变换精估计信号包括：
[0033]
s2-2-1、利用所述脑磁啁啾波变换粗估计信号与脑磁啁啾波变换原始数据计算粗估计结果的基础误差的计算式如下：
[0034][0035]
其中，e0为粗估计的基础误差，x为基础脑磁信号，a
emi
为第i次精估计对应的权重因子，g
emi
为第i次精估计对应的啁啾波，a
mpp
为粗估计得到的权重因子，g
mpp
为粗估计得到的啁啾波；
[0036]
s2-2-2、利用所述粗估计结果的基础误差计算脑磁啁啾波变换基础精细数据的计算式如下：
[0037][0038]
其中，yi为脑磁啁啾波变换基础精细数据，ai为残差信号权重因子，gi为第i次精细估计中的啁啾波，m为脑磁啁啾波变换数量，βi为误差系数因子e0为基础误差；
[0039]
s2-2-3、利用所述脑磁啁啾波变换基础精估计数据获取脑磁啁啾波变换基础精估计参数；
[0040]
s2-2-4、利用所述脑磁啁啾波变换基础精估计参数得到脑磁啁啾波变换精估计结果的残差权重因子与对应的高斯啁啾波；
[0041]
s2-2-5、利用所述脑磁啁啾变换波精估计的残差权重因子与高斯啁啾计算精细数据误差，所述精细数据误差的计算式如下：
[0042][0043]
其中，e为精估计误差，x为预处理后的原始脑磁信号，ai第i个脑磁啁啾波精细估计的残差权重因子，gi为第i个精估计的啁啾波；
[0044]
s2-2-6、判断所述精细数据误差是否小于精细数据阈值，若是，则利用所述精细数据误差对应的脑磁啁啾波变换基础精估计结果作为脑磁啁啾波变换精估计信号，否则，返回s2-2-1。
[0045]
进一步的，利用所述脑磁啁啾波变换精估计信号得到脑磁啁啾波变换信号包括：
[0046]
利用所述脑磁啁啾波变换精估计信号计算脑磁啁啾波变换精细数据残差数据的计算式如下：
[0047][0048]
其中，r
p
为精细估计后的脑磁啁啾波变换的残差，x为预处理后的原始脑磁信号，ai为残差权重因子，gi为第i个精估计的啁啾波；
[0049]
利用所述脑磁啁啾波变换精细数据残差数据计算脑磁啁啾波变换精估计相关系数的计算式如下：
[0050][0051]
其中，ccn为脑磁啁啾波变换精细数据相关系数，ai为残差权重因子，ri为第i个脑磁精细估计的残差，x为预处理后的原始脑磁信号，p-1为第p个啁啾波的索引；
[0052]
判断所述脑磁啁啾波变换精细数据残差数据对应的相关系数是否小于标准阈值，若是，则输出脑磁啁啾波变换精细数据残差数据对应的脑磁啁啾波变换精估计结果作为脑磁啁啾波变换基础数据，否则，返回s2-1；
[0053]
利用所述脑磁啁啾波变换基础数据计算脑磁啁啾波变换信号的计算式如下：
[0054][0055]
其中，为啁啾变换估计的脑磁信号，m为脑磁啁啾波变换数量，ai为残差权重因子，gi(t)为第i个脑磁啁啾波变换的啁啾波基；
[0056]
其中，标准阈值的范围为0至1。
[0057]
优选的，利用所述脑磁啁啾波变换信号得到基于啁啾波变换的脑磁啁啾波变换数据时频表征包括：
[0058]
利用所述脑磁啁啾波变换信号计算脑磁啁啾波变换wvd谱的计算式如下：
[0059][0060]
其中，t为时间，f为频率，gi为啁啾波，τ为脑磁信号的当地时间，g
i*
为啁啾波的共轭复数，e为自然底数，j为虚数的虚部；
[0061]
利用所述脑磁啁啾波变换wvd谱进行叠加处理得到脑磁啁啾波变换时频谱；
[0062]
利用所述脑磁啁啾波变换时频谱计算脑磁啁啾波变换数据时频表征的计算式如下：
[0063][0064]
其中，mpemact
x
为脑磁信号基于mp和em的啁啾变换的时频分布，t为时间，f为频率，m为脑磁啁啾波变换数量，ai为残差权重因子，gi啁啾变换得到的啁啾波；
[0065]
优选的，利用所述脑磁啁啾波变换数据时频表征得到脑磁信号瞬时激活频率包括：
[0066]
利用所述脑磁啁啾波变换数据时频表征基于局部最大值法得到脑磁信号瞬时激活频率。
[0067]
与最接近的现有技术相比，本发明具有的有益效果：
[0068]
通过精细线性调频波的数据估计，更加准确地对脑磁信号进行信号分离。同时，准确提供了有关脑磁信号的时间频率的能量分布信息，具有更高的时间频率分辨能力，能够同时较高分辨率地得到高频和低频时间分布特征。
附图说明
[0069]
图1是本发明提供的一种脑磁信号时频分布信息求解方法的流程图；
[0070]
图2是本发明提供的一种脑磁信号时频分布信息求解方法的实际应用整体流程图；
[0071]
图3是本发明提供的一种脑磁信号时频分布信息求解方法的实际应用线性调频波分解流程图；
[0072]
图4是本发明提供的一种脑磁信号时频分布信息求解方法的实际应用em算法流程图；
[0073]
图5是本发明提供的一种脑磁信号时频分布信息求解方法的效果对比图。
具体实施方式
[0074]
下面结合附图对本发明的具体实施方式作进一步的详细说明。
[0075]
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员
在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。
[0076]
实施例1：
[0077]
本发明提供了一种脑磁信号时频分布信息求解方法，如图1所示，包括：
[0078]
s1、获取初始脑磁信号进行降噪处理得到基础脑磁信号；
[0079]
s2、利用所述基础脑磁信号进行分解处理得到脑磁啁啾波变换信号；
[0080]
s3、利用所述脑磁啁啾变换波信号得到脑磁啁啾波变换数据时频表征；
[0081]
s4、利用所述脑磁啁啾波变换数据时频表征得到脑磁信号瞬时激活频率。
[0082]
s1具体包括：
[0083]
s1-1、获取初始脑磁信号进行降噪处理得到基础脑磁信号的计算式如下：
[0084][0085]
其中，x(t)为基础脑磁信号，t为时间，xk(t)是组成的第k个成分，n为成分个数，amk(t)为x(t)第k个成分的瞬时振幅，为x(t)第k个成分的瞬时相位，e为常数。
[0086]
s2具体包括：
[0087]
s2-1、利用所述基础脑磁信号基于mp算法得到脑磁啁啾波变换粗估计信号；
[0088]
s2-2、利用所述脑磁啁啾波变换粗估计信号基于em算法得到脑磁啁啾波变换精估计信号；
[0089]
s2-3、利用所述脑磁啁啾波变换精估计信号得到脑磁啁啾波变换信号。
[0090]
s2-1具体包括：
[0091]
s2-1-1、利用所述基础脑磁信号计算基础脑磁信号的残差信号权重因子的计算式如下：
[0092][0093]
其中，ai为残差信号权重因子，xi(t)为第i个基础脑磁信号，gi(t)为第i个啁啾变换的啁啾波，g
i*
(t)为gi(t)的共轭复数；
[0094]
s2-1-2、利用所述残差信号权重因子计算基础脑磁信号的残差信号的计算式如下：
[0095]
ri+1(t)＝xi(t)-aigi(t)
[0096]
其中，r
i+1
(t)为啁啾变换后的结果相对xi(t)的残差，xi(t)为第i个变换的基础脑磁信号成分，ai为残差信号权重因子，gi(t)为第i个啁啾变换的啁啾波；
[0097]
s2-1-3、利用所述残差信号计算残差信号能量的计算式如下：
[0098][0099]
其中，r
i+1
(t)为残差信号能量，r
i+1*
(t)为r
i+1
(t)的共轭复数；
[0100]
其中，用于脑磁啁啾波变换的啁啾波基的计算式如下：
[0101][0102]
其中，g(t)为脑磁啁啾波变换的啁啾波基，tc为啁啾波时间中心，fc为啁啾波的频
率中心，δt为采样间隔时间，c为啁啾率，d与δt相同，t为时间，j为虚数的虚部，i为虚数；
[0103]
s2-1-4、利用所述脑磁啁啾波变换的参数基于mp算法得到脑磁啁啾波变换粗估计信号；
[0104]
其中，脑磁啁啾波变换的估计参数为tc、fc、log(δt)、c、d。
[0105]
s2-2具体包括：
[0106]
s2-2-1、利用所述脑磁啁啾波变换粗估计信号与脑磁啁啾波变换原始数据计算粗估计结果的基础误差的计算式如下：
[0107][0108]
其中，e0为粗估计的基础误差，x为基础脑磁信号，a
emi
为第i次精估计对应的权重因子，g
emi
为第i次精估计对应的啁啾波，a
mpp
为粗估计得到的权重因子，g
mpp
为粗估计得到的啁啾波；
[0109]
s2-2-2、利用所述粗估计结果的基础误差计算脑磁啁啾波变换基础精细数据的计算式如下：
[0110][0111]
其中，yi为脑磁啁啾波变换基础精细数据，ai为残差信号权重因子，gi为第i次精细估计中的啁啾波，m为脑磁啁啾波变换数量，βi为误差系数因子e0为基础误差；
[0112]
s2-2-3、利用所述脑磁啁啾波变换基础精估计数据获取脑磁啁啾波变换基础精估计参数；
[0113]
s2-2-4、利用所述脑磁啁啾波变换基础精估计参数得到脑磁啁啾波变换精估计结果的残差权重因子与对应的高斯啁啾波；
[0114]
s2-2-5、利用所述脑磁啁啾变换波精估计的残差权重因子与高斯啁啾计算精细数据误差，所述精细数据误差的计算式如下：
[0115][0116]
其中，e为精估计误差，x为预处理后的原始脑磁信号，ai第i个脑磁啁啾波精细估计的残差权重因子，gi为第i个精估计的啁啾波；
[0117]
s2-2-6、判断所述精细数据误差是否小于精细数据阈值，若是，则利用所述精细数据误差对应的脑磁啁啾波变换基础精估计结果作为脑磁啁啾波变换精估计信号，否则，返回s2-2-1。
[0118]
s2-3具体包括：
[0119]
s2-3-1、利用所述脑磁啁啾波变换精估计信号计算脑磁啁啾波变换精细数据残差数据的计算式如下：
[0120][0121]
其中，r
p
为精细估计后的脑磁啁啾波变换的残差，x为预处理后的原始脑磁信号，ai为残差权重因子，gi为第i个精估计的啁啾波；
[0122]
s2-3-2、利用所述脑磁啁啾波变换精细数据残差数据计算脑磁啁啾波变换精估计相关系数的计算式如下：
[0123][0124]
其中，ccn为脑磁啁啾波变换精细数据相关系数，ai为残差权重因子，ri为第i个脑磁精细估计的残差，x为预处理后的原始脑磁信号，p-1为第p个啁啾波的索引；
[0125]
s2-3-3、判断所述脑磁啁啾波变换精细数据残差数据对应的相关系数是否小于标准阈值，若是，则输出脑磁啁啾波变换精细数据残差数据对应的脑磁啁啾波变换精估计结果作为脑磁啁啾波变换基础数据，否则，返回s2-1；
[0126]
s2-3-4、利用所述脑磁啁啾波变换基础数据计算脑磁啁啾波变换信号的计算式如下：
[0127][0128]
其中，为啁啾变换估计的脑磁信号，m为脑磁啁啾波变换数量，ai为残差权重因子，gi(t)为第i个脑磁啁啾波变换的啁啾波基；
[0129]
其中，标准阈值的范围为0至1。
[0130]
s3具体包括：
[0131]
s3-1、利用所述脑磁啁啾波变换信号计算脑磁啁啾波变换wvd谱的计算式如下：
[0132][0133]
其中，t为时间，f为频率，gi为啁啾波，τ为脑磁信号的当地时间，g
i*
为啁啾波的共轭复数，e为自然底数，j为虚数的虚部；
[0134]
s3-2、利用所述脑磁啁啾波变换wvd谱进行叠加处理得到脑磁啁啾波变换时频谱；
[0135]
s3-3、利用所述脑磁啁啾波变换时频谱计算脑磁啁啾波变换数据时频表征的计算式如下：
[0136][0137]
其中，mpemact
x
为脑磁信号基于mp和em的啁啾变换的时频分布，t为时间，f为频率，m为脑磁啁啾波变换数量，ai为残差权重因子，gi啁啾变换得到的啁啾波；
[0138]
s4具体包括：
[0139]
s4-1、利用所述脑磁啁啾波变换数据时频表征基于局部最大值法得到脑磁信号瞬
时激活频率。
[0140]
实施例2：
[0141]
本发明提供了一种脑磁信号时频分布信息求解方法的实际应用流程，如图2所示，包括：
[0142]
通过用稀疏且能量保持的高斯啁啾信号来表征脑磁信号，同时利用该表征可以求得脑磁信号的时频分布，以及瞬时频率激活情况。
[0143]
步骤11为利用基于原子磁强计阵列测得，并经过预处理降噪的脑磁信号x(t)，具有多成分分量，定义如下：
[0144][0145]
其中，t为时间，xk(t)是组成x(t)的第k个成分，n为成分个数，其中amk(t)是x(t)第k个成分的瞬时振幅；是x(t)第k个成分的瞬时相位，e为常数。
[0146]
步骤12为基于匹配追踪(mp)和期望最大值(em)对脑磁信号进行chirplet分解，将脑磁信号分解为多个chirplet波的总和。
[0147]
chirplet波选择的是高斯啁啾，其一般形式表达式为
[0148][0149]
本方法中需要优化估计的参数为i＝{tc，fc，log(δt)，c，d}，tc为chirplet波时间中心，fc为chirplet波的频率中心，δt为采样间隔时间，c为chirplet率，d＝δt，t为时间。
[0150]
x(t)基于g(t)的chirplet变换如下：
[0151][0152]
其中，g
*
(t)为g(t)的共轭。
[0153]
步骤13为对步骤12中求得的各个chirplert波求wvd谱的计算式如下：
[0154][0155]
步骤14为对步骤13中求得的wvd谱进行叠加，得到chirplet变换时频谱，从而得到chirplet的时频表征的计算式如下：
[0156][0157]
步骤15为利用局部最大值的方法对步骤14求得的时频谱进行计算，最后求得步骤16瞬时激活频谱。
[0158]
步骤12chirplet分解的具体实现方式，如图3所示，具体实现方式为：
[0159]
步骤21为根据表1初始化第一个chirplet波初始阶数p＝0，初始化信号残差r0＝x(t)；
[0160]
步骤22为利用mp算法粗估计第i个chirplet，即找到参数ii＝{tc，fc，log(δt)，c，
d}的最优解，并自适应地生成gi(t)以最小化剩余信号的能量，具体实现方法为：
[0161]
残差信号为
[0162]ri+1
(t)＝xi(t)-a
igi
(t) (25)
[0163]
初始值为：i＝0，r0＝x(t)。
[0164]
其中，权重因子ai为：
[0165][0166]
残差信号能量为：
[0167][0168]
当||r
i+1
(t)||2最小时，可以得到第p个chirplet的mp粗估计的参数为i
mpp
＝{tc，fc，log(δt)，c，d}
mp
，对应的chirplet和权重因子分别为g
mpp
和a
mpp
。
[0169]
步骤23为利用em方法估计对mp粗估计后的第i个chirplet进行精估计。
[0170]
步骤24的具体实现方式为：由步骤23得到精确估计的gi和ai后，重新计算残差
[0171][0172]
残差若满足的条件——相关系数ccn小于设定值(该值取值范围为[0,1])时，则停止chirplet分解，若不满足重复步骤22-24。
[0173][0174]
步骤25为得到的由采集的脑磁信号x(t)估计的m个chirplert波，则有脑磁信号的chirplet变换为
[0175][0176]
其中步骤23中的em算法的具体实现步骤，如图4所示：步骤31为e-step，计算原始数据和已估计结果之间的误差：
[0177][0178]
可以得到需要进一步精估计的数据为：
[0179][0180]
其中，i＝1,2,3,
…
,p
[0181][0182]
步骤32为m-step，是利用mp对yi进行重新估计参数i
emi
＝{tc，fc，log(δt)，c，d)
em
，
更新gi和ai。
[0183]
步骤33为重新计算原始数据和估计结果之间的误差：
[0184][0185]
重复步骤31-33，当e小于设定阈值时，就能得到精确估计的gi和ai。
[0186]
上述方案中提到的各值定义如下表：
[0187][0188]
[0189]
本实施例中，一种脑磁信号时频分布信息求解方法的实际应用流程，如图5所示，该自适应的chirplet变换技术中，通过22-24步骤的精细chirplet估计，更加准确地对脑磁信号进行信号分离。同时利用13-16步骤，准确提供了有关脑磁信号的时间频率的能量分布信息，相较于stft、cwt具有更高的时间频率分辨能力，能够同时较高分辨率地得到高频和低频时间分布特征。
[0190]
最后应当说明的是：以上实施例仅用以说明本发明的技术方案而非对其限制，尽管参照上述实施例对本发明进行了详细的说明，所属领域的普通技术人员应当理解：依然可以对本发明的具体实施方式进行修改或者等同替换，而未脱离本发明精神和范围的任何修改或者等同替换，其均应涵盖在本发明的权利要求保护范围之内。