CN113951887A - 一种频谱匹配独立成分分析方法及系统 - Google Patents

一种频谱匹配独立成分分析方法及系统 Download PDF

Info

Publication number
CN113951887A
CN113951887A CN202111583666.0A CN202111583666A CN113951887A CN 113951887 A CN113951887 A CN 113951887A CN 202111583666 A CN202111583666 A CN 202111583666A CN 113951887 A CN113951887 A CN 113951887A
Authority
CN
China
Prior art keywords
subject data
meg system
meg
acquired
frequency band
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202111583666.0A
Other languages
English (en)
Other versions
CN113951887B (zh
Inventor
高阳
李勇
宁晓琳
房建成
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hangzhou Innovation Research Institute of Beihang University
Original Assignee
Hangzhou Innovation Research Institute of Beihang University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hangzhou Innovation Research Institute of Beihang University filed Critical Hangzhou Innovation Research Institute of Beihang University
Priority to CN202111583666.0A priority Critical patent/CN113951887B/zh
Publication of CN113951887A publication Critical patent/CN113951887A/zh
Application granted granted Critical
Publication of CN113951887B publication Critical patent/CN113951887B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/24—Detecting, measuring or recording bioelectric or biomagnetic signals of the body or parts thereof
    • A61B5/242—Detecting biomagnetic fields, e.g. magnetic fields produced by bioelectric currents
    • A61B5/245—Detecting biomagnetic fields, e.g. magnetic fields produced by bioelectric currents specially adapted for magnetoencephalographic [MEG] signals
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7203—Signal processing specially adapted for physiological signals or for diagnostic purposes for noise prevention, reduction or removal
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7225—Details of analog processing, e.g. isolation amplifier, gain or sensitivity adjustment, filtering, baseline or drift compensation
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7235—Details of waveform analysis
    • A61B5/7253—Details of waveform analysis characterised by using transforms
    • A61B5/7257—Details of waveform analysis characterised by using transforms using Fourier transforms
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
    • G06F18/2134—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on separation criteria, e.g. independent component analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/22—Matching criteria, e.g. proximity measures

Landscapes

  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Artificial Intelligence (AREA)
  • Signal Processing (AREA)
  • Medical Informatics (AREA)
  • Surgery (AREA)
  • Veterinary Medicine (AREA)
  • Public Health (AREA)
  • General Health & Medical Sciences (AREA)
  • Theoretical Computer Science (AREA)
  • Animal Behavior & Ethology (AREA)
  • Molecular Biology (AREA)
  • Data Mining & Analysis (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Biophysics (AREA)
  • Pathology (AREA)
  • Biomedical Technology (AREA)
  • Psychiatry (AREA)
  • Physiology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Evolutionary Computation (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Biology (AREA)
  • Power Engineering (AREA)
  • Mathematical Physics (AREA)
  • Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)

Abstract

本发明涉及一种频谱匹配独立成分分析方法及系统，方法包括：获取MEG系统所采集到的被试者数据中单个频带内的噪声功率B a 参数的具体值、单个频带内的源功率A a 参数的具体值、被试者数据频域内的噪声项的频谱协方差矩阵和被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值；基于MEG系统所采集到的被试者数据中单个频带内的噪声功率B a 参数的具体值、单个频带内的源功率A a 参数的具体值、被试者数据频域内的噪声项的频谱协方差矩阵和被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值，获取该MEG系统所采集到的被试者数据中预先设定的频带内的进行滤波的信号矩阵X a 中的源信号。

Description

一种频谱匹配独立成分分析方法及系统
技术领域
本发明涉及脑磁图医疗设备数据处理技术领域，尤其涉及一种频谱匹配独立成分分析方法及系统。
背景技术
由于脑磁图具有高时间分辨率与空间分辨率特性，其在脑科学研究和大脑疾病诊断等方面具有广泛应用前景。然而，现有的基于超导量子干涉装置(SQUID)的脑磁图系统需要在低温下运行，因此必须固定在低温杜瓦瓶内，使得系统变得庞大、笨重且它们无法适应不同的头部形状或尺寸，为脑磁图系统的信号测量精度和使用场景带来了巨大的挑战。近年来量子技术的进步促进了新型磁场传感器的发展，基于量子传感技术的光泵磁强计(OPM)可提供亚fT级极弱磁磁场测量，其灵敏度与传统MEG中使用的SQUID传感器相似，但它们不需要冷却，可贴近头皮进行测量，而且可以集成芯片化,其体积轻便，因此基于OPM（Optical Pump Magnetometer：光泵磁强计）的脑磁图系统被广泛的研究。
目前基于OPM的可穿戴MEG（Magnetoencephalogram：脑磁图）系统已经被研制出来，其具备能够提供较高的数据质量、适应任何受试者、允许在移动过程中捕获数据等特点，已经被许多研究机构引入到脑科学研究领域进行信号测量。但是，应用基于OPM的可穿戴MEG系统进行信号测量分析，关键的步骤是进行源分离，因为测到的信号中不仅包含了大脑活动信号源，还会包含眨眼、心跳、线噪声、肌肉活动等干扰源，这些干扰源的存在会对源定位和分析产生较大影响，导致分析结果不准确甚至出错。
在脑磁图信号分析领域进行源分离较为流行的方法是独立成分分析，但现存的独立成分分析方法大多采用非高斯独立源的无噪声混合模型，其存在强制估计与传感器一样多的源，依赖于损失信号的PCA降维操作等限制条件，导致其对噪声变化较大的移动情况下捕捉到的数据进行源分离准确性较差，使得基于OPM的脑磁图不能较好的应用于人体移动情况下的数据测量与分析。
发明内容
（一）要解决的技术问题
鉴于现有技术的上述缺点、不足，本发明提供一种频谱匹配独立成分分析方法及系统，其解决了现有的独立成分分析方法大多采用非高斯独立源的无噪声混合模型，其存在强制估计与传感器一样多的源，依赖于损失信号的PCA降维操作等限制条件，导致其对噪声变化较大的移动情况下捕捉到的数据进行源分离准确性较差技术问题。
（二）技术方案
为了达到上述目的，本发明采用的主要技术方案包括：
第一方面，本发明实施例提供一种频谱匹配独立成分分析方法，包括：
S1、获取MEG系统所采集到的被试者数据中单个频带内的噪声功率B a 参数的具体值、单个频带内的源功率A a 参数的具体值、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值；
所述噪声功率B a 参数为MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵在预先划分的频带中的第a个频带内的对应平均值；
所述源功率A a 参数为MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵在预先划分的频带中的第a个频带内的对应平均值；
S2、基于所述MEG系统所采集到的被试者数据中单个频带内的噪声功率B a 参数的具体值、单个频带内的源功率A a 参数的具体值、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值，获取该MEG系统所采集到的被试者数据中预先设定的频带内的进行滤波的信号矩阵X a 中的源信号。
优选的，所述S2具体包括：
基于所述MEG系统所采集到的被试者数据中单个频带内的噪声功率B a 参数的具体值、单个频带内的源功率A a 参数的具体值、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值，根据公式（1）计算得到该MEG系统所采集到的被试者数据中第a频带范围上的进行滤波的信号矩阵X a 中的源信号；
所述公式（1）为：
Figure DEST_PATH_IMAGE001
其中，
Figure 527264DEST_PATH_IMAGE002
为该MEG系统所采集到的被试者数据中预先划分的频带中的第a个频带内 的进行滤波的信号矩阵X a 对应恢复的源信号矩阵；
Figure DEST_PATH_IMAGE003
为MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系 统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的转置；
Figure 523033DEST_PATH_IMAGE004
为所述MEG系统所采集到的被试者数据中单个频带内的噪声功率B a 参数的逆；
Figure DEST_PATH_IMAGE005
为单个频带内的源功率A a 参数的逆。
优选的，在S1之前还包括：
S0、基于预先获取的MEG系统所采集到的被试者数据、预先获取的背景噪音信号，建立频谱模型，所述频谱模型中包括噪声功率B a 参数、单个频带内的源功率A a 参数、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M参数；
相应的，所述S1具体包括：
基于预先获取的MEG系统所采集到的被试者数据、预先获取的MEG系统采集的空房间中的背景噪声信号和所述频谱模型，获取MEG系统所采集到的被试者数据中单个频带内的噪声功率B a 参数的具体值、单个频带内的源功率A a 参数的具体值、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值。
优选的，
所述频谱模型为：
Figure 630666DEST_PATH_IMAGE006
其中，C a 为C（f）在所划分的单个频带内的对应平均值；
C（f）为X（t）的自协方差矩阵的傅里叶变换；
其中，所述X（t）为p维向量，表示从具有p个探头传感器的MEG系统所采集到的被试者数据；
其中，
Figure DEST_PATH_IMAGE007
其中，N（t）为MEG系统所采集到的被试者数据时域内的背景噪声项且为p维向量；
Figure 337460DEST_PATH_IMAGE008
表示MEG系统所采集到的被试者数据时域内的q个独立源且为q维向量；
M是一个
Figure DEST_PATH_IMAGE009
的混合矩阵，表示假设通过该矩阵源信号与噪声信号实现加性混合。
优选的，所述S1具体包括：
S11、基于预先获取的MEG系统所采集到的被试者数据，获取C a 参数的具体值；
S12、基于预先获取的背景噪声信号，获取B a 参数的具体值；
S13、基于所述C a 参数的具体值和所述B a 参数的具体值，获取单个频带内的源功率A a 参数的具体值、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值。
优选的，所述S11具体包括：
基于预先获取的MEG系统所采集到的被试者数据，采用公式（2）获取预先获取的MEG系统所采集到的被试者数据X（t）的离散样本X（r）的傅里叶系数；
所述公式（2）为：
Figure 214149DEST_PATH_IMAGE010
；
Figure DEST_PATH_IMAGE011
为MEG系统所采集到的被试者数据X（t）的离散样本X（r）的傅里叶系数；
R表示预先获取的MEG系统所采集到的被试者数据时域采样点数；r表示离散样本X （r）中的第r个样本；
基于所述X（r）的傅里叶系数
Figure 849661DEST_PATH_IMAGE011
，采用公式（3）获取
Figure 647852DEST_PATH_IMAGE012
的具体值，并将所述
Figure 377911DEST_PATH_IMAGE012
的具 体值作为C a 的具体值；
所述公式（3）为：
Figure DEST_PATH_IMAGE013
；
Figure 904181DEST_PATH_IMAGE012
为C a 的估计值；
H表示对矩阵做共轭转置；
Figure 694283DEST_PATH_IMAGE014
为
Figure 730503DEST_PATH_IMAGE011
的共轭转置；
其中，n a 为预先划分MEG系统所采集到的被试者数据的频带中的第a个频带内可用的傅里叶系数数目，并且，
Figure DEST_PATH_IMAGE015
；
其中，#表示计数；
I a 为预先划分的频带中的第a个频带；
k表示总的傅里叶系数数目。
优选的，所述S12具体包括：
基于预先获取的MEG系统所采集到空房间中的背景噪声信号，采用公式（4）获取预先获取的MEG系统所采集到的背景噪声信号N（t）的离散样本N（w）的傅里叶系数；
所述公式（4）为：
Figure 529832DEST_PATH_IMAGE016
；
Figure DEST_PATH_IMAGE017
为MEG系统所采集到的背景噪声信号N（t）的离散样本N（w）的傅里叶系数；
W表示MEG系统所采集到的背景噪声信号时域采样点数；
w表示离散样本N（w）中的第w个样本；
基于所述N（w）的傅里叶系数
Figure 318796DEST_PATH_IMAGE017
，采用公式（5）获取
Figure 529066DEST_PATH_IMAGE018
的具体值，并将所述
Figure 567430DEST_PATH_IMAGE018
的具 体值作为B a 的具体值；
所述公式（5）为：
Figure DEST_PATH_IMAGE019
；
Figure 904870DEST_PATH_IMAGE018
为B a 的估计值；
Figure 282762DEST_PATH_IMAGE020
为
Figure 165398DEST_PATH_IMAGE017
的共轭转置；
其中，z a 为预先划分MEG系统所采集到的背景噪声信号的频带中的第a个频带内可用的傅里叶系数数目，并且，
Figure DEST_PATH_IMAGE021
。
优选的，所述S13具体包括：
基于KL散度构造频谱匹配标准，所述频谱匹配标准的形式为：
Figure 159899DEST_PATH_IMAGE022
其中，m为频带数目；
Figure DEST_PATH_IMAGE023
表示
Figure 566610DEST_PATH_IMAGE024
和
Figure DEST_PATH_IMAGE025
的KL散度；
基于所述频谱匹配标准、所述C a 参数的具体值和所述B a 参数的具体值，采用期望最大化算法获取单个频带内的源功率A a 参数的具体值和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值。
优选的，
所述MEG系统所采集到的被试者数据为用户处于背景剩磁消减后的磁屏蔽房内并在相同声音信号的刺激下静止坐下、静止站立以及身体和头部移动时所采集的MEG被试者数据；
所述背景剩磁消减后的磁屏蔽房内的背景磁场低于50nT；
预先获取的背景噪音信号为背景剩磁消减后的磁屏蔽房内的预先设定时间段内的空的磁屏蔽房内的磁信号。
第二方面，本发明实施例提供一种频谱匹配独立成分分析系统，所述系统包括：
至少一个处理器；以及
与所述处理器通信连接的至少一个存储器，其中，所述存储器存储有可被所述处理器执行的程序指令，所述处理器调用所述程序指令能够执行上述任一的频谱匹配独立成分分析方法。
（三）有益效果
本发明的有益效果是：本发明的一种频谱匹配独立成分分析方法及系统，由于通过获取的MEG系统所采集到的被试者数据中单个频带内的噪声功率B a 参数的具体值、单个频带内的源功率参数A a 的具体值、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值，得到该MEG系统所采集到的被试者数据中预先设定的频带内的进行滤波的信号矩阵X a 中的源信号，使得源信号恢复更为准确。
附图说明
图1为本发明的一种频谱匹配独立成分分析方法流程图；
图2为本发明实施例二中的MEG系统中的磁屏蔽系统结构示意图；
图3为本发明实施例二中的MEG系统中的OPM-MEG头盔阵列数据获取系统结构示意图；
图4为本发明实施例二中的MEG系统中的听力刺激系统结构示意图；
图5为本发明实施例二中的MEG系统在采集被试者数据时的示意图；
图6为本发明实施例二中的频谱匹配独立成分分析方法流程图。
【附图标记说明】
101：磁屏蔽房；
102：磁场消零线圈；
103：参考阵列；
104：线圈驱动模块；
105：NI DAQ；
1051：DAC；
1052：ADC；
106：传感器采集控制电路；
107：传感器、线圈控制设备；
201：头盔阵列；
203：数据获取设备；
301：无磁性塑料细管耳机；
302：听力刺激设备。
具体实施方式
为了更好的解释本发明，以便于理解，下面结合附图，通过具体实施方式，对本发明作详细描述。
为了更好的理解上述技术方案，下面将参照附图更详细地描述本发明的示例性实施例。虽然附图中显示了本发明的示例性实施例，然而应当理解，可以以各种形式实现本发明而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了能够更清楚、透彻地理解本发明，并且能够将本发明的范围完整的传达给本领域的技术人员。
实施例一
参见图1，本实施例提供一种频谱匹配独立成分分析方法，包括：
S1、获取MEG系统所采集到的被试者数据中单个频带内的噪声功率B a 参数的具体值、单个频带内的源功率A a 参数的具体值、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值。
所述噪声功率B a 参数为MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵在预先划分的频带中的第a个频带内的对应平均值。
所述源功率A a 参数为MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵在预先划分的频带中的第a个频带内的对应平均值。
具体的，本实施例中MEG系统所采集到的被试者数据的频率范围为0-70HZ，然后被预先划分为40个频带，每个频带的频率间隔范围都相同。
S2、基于所述MEG系统所采集到的被试者数据中单个频带内的噪声功率B a 参数的具体值、单个频带内的源功率A a 参数的具体值、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值，获取该MEG系统所采集到的被试者数据中预先设定的频带内的进行滤波的信号矩阵X a 中的源信号。
在本实施例的实际应用中，所述S2具体包括：
基于所述MEG系统所采集到的被试者数据中单个频带内的噪声功率B a 参数的具体值、单个频带内的源功率A a 参数的具体值、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值，根据公式（1）计算得到该MEG系统所采集到的被试者数据中第a频带范围上的进行滤波的信号矩阵X a 中的源信号。
所述公式（1）为：
Figure 582364DEST_PATH_IMAGE001
其中，
Figure 885169DEST_PATH_IMAGE002
为该MEG系统所采集到的被试者数据中预先划分的频带中的第a个频带内 的进行滤波的信号矩阵X a 对应恢复的源信号矩阵。
Figure 632545DEST_PATH_IMAGE003
为MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系 统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的转置。
Figure 46209DEST_PATH_IMAGE004
为所述MEG系统所采集到的被试者数据中单个频带内的噪声功率B a 参数的逆。
Figure 680584DEST_PATH_IMAGE005
为单个频带内的源功率A a 参数的逆。
在本实施例的实际应用中，在S1之前还包括：
S0、基于预先获取的MEG系统所采集到的被试者数据、预先获取的背景噪音信号，建立频谱模型，所述频谱模型中包括噪声功率B a 参数、单个频带内的源功率A a 参数、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M参数。
相应的，所述S1具体包括：
基于预先获取的MEG系统所采集到的被试者数据、预先获取的MEG系统采集的空房间中的背景噪声信号和所述频谱模型，获取MEG系统所采集到的被试者数据中单个频带内的噪声功率B a 参数的具体值、单个频带内的源功率A a 参数的具体值、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值。
在本实施例的实际应用中，所述频谱模型为：
Figure 357553DEST_PATH_IMAGE006
其中，C a 为C（f） 在所划分的单个频带内的对应平均值；C（f）为X（t）的自协方差矩阵的傅里叶变换；其中，所 述X（t）为p维向量，表示从具有p个探头传感器的MEG系统所采集到的被试者数据。
其中，
Figure 326646DEST_PATH_IMAGE007
其中，N（t）为MEG系统所采集到的被试者数据时域内的 背景噪声项且为p维向量；
Figure 544000DEST_PATH_IMAGE008
表示MEG系统所采集到的被试者数据时域内的q个独立源且 为q维向量；M是一个
Figure 547728DEST_PATH_IMAGE009
的混合矩阵，表示假设通过该矩阵源信号与噪声信号实现加性混 合。
在本实施例的实际应用中，所述S1具体包括：
S11、基于预先获取的MEG系统所采集到的被试者数据，获取C a 参数的具体值。
S12、基于预先获取的背景噪声信号，获取B a 参数的具体值。
S13、基于所述C a 参数的具体值和所述B a 参数的具体值，获取单个频带内的源功率A a 参数的具体值、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值。
在本实施例的实际应用中，所述S11具体包括：
基于预先获取的MEG系统所采集到的被试者数据，采用公式（2）获取预先获取的MEG系统所采集到的被试者数据X（t）的离散样本X（r）的傅里叶系数。
所述公式（2）为：
Figure 395599DEST_PATH_IMAGE010
；
Figure 897993DEST_PATH_IMAGE011
为MEG系统所采集到的被试者数据X（t）的离散样本X（r）的傅里叶系数。R表示 预先获取的MEG系统所采集到的被试者数据时域采样点数；r表示离散样本X（r）中的第r个 样本。
基于所述X（r）的傅里叶系数
Figure 856722DEST_PATH_IMAGE011
，采用公式（3）获取
Figure 449377DEST_PATH_IMAGE012
的具体值，并将所述
Figure 264886DEST_PATH_IMAGE012
的具 体值作为C a 的具体值。
所述公式（3）为：
Figure 959304DEST_PATH_IMAGE013
。
Figure 518461DEST_PATH_IMAGE012
为C a 的估计值。
H表示对矩阵做共轭转置；
Figure 231203DEST_PATH_IMAGE014
为
Figure 217613DEST_PATH_IMAGE011
的共轭转置。
其中，n a 为预先划分MEG系统所采集到的被试者数据的频带中的第a个频带内可用的傅里叶系数数目，并且，
Figure 963109DEST_PATH_IMAGE015
；
其中，#表示计数。
I a 为预先划分的频带中的第a个频带。
k表示总的傅里叶系数数目。
在本实施例的实际应用中，所述S12具体包括：
基于预先获取的MEG系统所采集到空房间中的背景噪声信号，采用公式（4）获取预先获取的MEG系统所采集到的背景噪声信号N（t）的离散样本N（w）的傅里叶系数。
所述公式（4）为：
Figure 60378DEST_PATH_IMAGE016
；
Figure 643937DEST_PATH_IMAGE017
为MEG系统所采集到的背景噪声信号N（t）的离散样本N（w）的傅里叶系数。
W表示MEG系统所采集到的背景噪声信号时域采样点数。
w表示离散样本N（w）中的第w个样本。
基于所述N（w）的傅里叶系数
Figure 738932DEST_PATH_IMAGE017
，采用公式（5）获取
Figure 922788DEST_PATH_IMAGE018
的具体值，并将所述
Figure 823748DEST_PATH_IMAGE018
的具 体值作为B a 的具体值。
所述公式（5）为：
Figure 245502DEST_PATH_IMAGE019
。
Figure 776978DEST_PATH_IMAGE018
为B a 的估计值。
Figure 166240DEST_PATH_IMAGE020
为
Figure 870891DEST_PATH_IMAGE017
的共轭转置。
其中，z a 为预先划分MEG系统所采集到的背景噪声信号的频带中的第a个频带内可用的傅里叶系数数目，并且，
Figure 943889DEST_PATH_IMAGE021
。
在本实施例的实际应用中，所述S13具体包括：
基于KL散度构造频谱匹配标准，所述频谱匹配标准的形式为：
Figure 459315DEST_PATH_IMAGE022
其中，m为频带数目。
Figure 289868DEST_PATH_IMAGE023
表示
Figure 532630DEST_PATH_IMAGE024
和
Figure 194556DEST_PATH_IMAGE025
的KL散度。
基于所述频谱匹配标准、所述C a 参数的具体值和所述B a 参数的具体值，采用期望最大化算法获取单个频带内的源功率A a 参数的具体值和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值。
在本实施例的实际应用中，所述MEG系统所采集到的被试者数据为用户处于背景剩磁消减后的磁屏蔽房内并在相同声音信号的刺激下静止坐下、静止站立以及身体和头部移动时所采集的MEG被试者数据。
所述背景剩磁消减后的磁屏蔽房内的背景磁场低于50nT。
预先获取的背景噪音信号为背景剩磁消减后的磁屏蔽房内的预先设定时间段内的空的磁屏蔽房内的磁信号。
另一方面，本实施例还提供一种频谱匹配独立成分分析系统，所述系统包括：
至少一个处理器；以及与所述处理器通信连接的至少一个存储器，其中，所述存储器存储有可被所述处理器执行的程序指令，所述处理器调用所述程序指令能够执行上述频谱匹配独立成分分析方法。
由于本发明上述实施例所描述的系统，为实施本发明上述实施例的方法所采用的系统，故而基于本发明上述实施例所描述的方法，本领域所属技术人员能够了解该系统的具体结构及变形，因而在此不再赘述。凡是本发明上述实施例的方法所采用的系统都属于本发明所欲保护的范围。
本实施例中的一种频谱匹配独立成分分析方法及系统，由于通过获取的MEG系统所采集到的被试者数据中单个频带内的噪声功率B a 参数的具体值、单个频带内的源功率参数A a 的具体值、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值，得到该MEG系统所采集到的被试者数据中预先设定的频带内的进行滤波的信号矩阵X a 中的源信号，使得源信号恢复更为准确。
实施例二
参见图2-5，本实施例中的MEG系统包括：磁屏蔽系统、OPM-MEG头盔阵列数据获取系统、听力刺激系统。
其中，本实施例中磁屏蔽系统包括磁屏蔽房101、磁场消零线圈102、由XYZ方向构成的参考阵列103、NI DAQ（NI数据采集系统）105、传感器采集控制电路106、传感器和线圈控制设备107、线圈驱动模块104等构件。
其中，NI DAQ105包括：DAC（数模转换器）1051，ADC（模数转换器）1052。
OPM-MEG头盔阵列数据获取系统包括头盔阵列201（由3D打印头盔和一套32通道OPM探头组成）、NI DAQ105、数据获取设备203等构件。
听力刺激系统包括听力刺激设备302、无磁性塑料细管耳机301等构件。
MEG系统采集信号前，会启动磁屏蔽房进行主动消磁，在完成主动消磁后打开OPM传感器进行参考信号测量，利用测量的参考信号和剩磁消零算法进行背景剩磁消减，以实现最大程度的背景磁场抑制。
开始进行空房间的磁信号测量，记录时长与后面被试者进入房间到完成实验的计划时长一致（也就是预先设定时间段），将其视为背景噪声信号。
然后 ，受试者进入房间，背上提前准备好的双肩背包，戴上3D打印的无磁性头盔，并将头盔线固定在受试者所背的双肩背包上用夹子夹上。根据实验需要将32个探头分别插入对应插槽，上下两层排列阵列完全相同，最后将塑料细管耳机戴在耳朵上。
最后关上磁屏蔽房门，再重复上述剩磁补偿步骤，至此所述背景剩磁消减后的磁屏蔽房内的背景磁场低于50nT，完成数据采集前所有准备。
发出声音信号，指示受试者按声音信号要求做出相应动作，并播放声音信号，声音信号由高保真声卡通过无磁性塑料细管耳机呈现给受试者。同样的声音信号的刺激作用下，受试者身体状况在静止坐下、静止站立、身体和头部移动情况下进行数据采集，得到MEG系统所采集到的被试者数据。
本实施例中通过测量空房间下的信号作为背景噪声信号，使得模型解算更为接近实际情况，源分离准确度提高。
本实施例在声音刺激下采用了静态坐立、静态站立、身体和大脑移动等三种状态进行数据采集，增强了数据的可参照性。
本实施例中MEG系统的数据采集能够实现较高质量自然状态的数据采集，使得测量数据具有可对比性和参照性。
参见图6，本实施例中的频谱匹配独立成分分析方法包括：
Q0、基于预先获取的MEG系统所采集到的被试者数据、预先获取的背景噪音信号，建立频谱模型，所述频谱模型中包括噪声功率B a 参数、单个频带内的源功率A a 参数、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M参数。
所述频谱模型为：
Figure 864571DEST_PATH_IMAGE006
其中，C a 为C（f）在所划分的单个频带内的对应平均值。
C（f）为X（t）的自协方差矩阵的傅里叶变换。
其中，所述X（t）为p维向量，表示从具有p个探头传感器的MEG系统所采集到的被试者数据。
其中，
Figure 182420DEST_PATH_IMAGE007
其中，N（t）为MEG系统所采集到的被试者数据时域内的背景噪声项且为p维向量。
Figure 215492DEST_PATH_IMAGE008
表示MEG系统所采集到的被试者数据时域内的q个独立源且为q维向量。
M是一个
Figure 997503DEST_PATH_IMAGE009
的混合矩阵，表示假设通过该矩阵源信号与噪声信号实现加性混合。
本实施例中的频谱模型能够将预先获取的MEG系统所采集到的被试者数据降维嵌套入模型，克服了传统独立成分分析需要假设源与噪声相同可能导致的模型退化问题。频谱模型构造假设源和噪声信号为平稳的高斯时间序列，考虑到了源和噪声信号的时间相关性，使得源信号恢复更准确。
Q1、基于预先获取的MEG系统所采集到的被试者数据、预先获取的MEG系统采集的空房间中的背景噪声信号和所述频谱模型，获取MEG系统所采集到的被试者数据中单个频带内的噪声功率B a 参数的具体值、单个频带内的源功率A a 参数的具体值、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值。
所述噪声功率B a 参数为MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵在预先划分的频带中的第a个频带内的对应平均值。
所述源功率A a 参数为MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵在预先划分的频带中的第a个频带内的对应平均值。
Q2、基于所述MEG系统所采集到的被试者数据中单个频带内的噪声功率B a 参数的具体值、单个频带内的源功率A a 参数的具体值、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值，根据公式（1）计算得到该MEG系统所采集到的被试者数据中第a频带范围上的进行滤波的信号矩阵X a 中的源信号。
所述公式（1）为：
Figure 838420DEST_PATH_IMAGE001
其中，
Figure 643565DEST_PATH_IMAGE002
为该MEG系统所采集到的被试者数据中预先划分的频带中的第a个频带内 的进行滤波的信号矩阵X a 对应恢复的源信号矩阵。
Figure 228130DEST_PATH_IMAGE003
为MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系 统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的转置。
Figure 615380DEST_PATH_IMAGE004
为所述MEG系统所采集到的被试者数据中单个频带内的噪声功率B a 参数的逆。
Figure 627199DEST_PATH_IMAGE005
为单个频带内的源功率A a 参数的逆。
具体的，在本实施例中，所述Q1具体包括：
Q11、基于预先获取的MEG系统所采集到的被试者数据，获取C a 参数的具体值。
Q12、基于预先获取的背景噪声信号，获取B a 参数的具体值。
Q13、基于所述C a 参数的具体值和所述B a 参数的具体值，获取单个频带内的源功率A a 参数的具体值、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值。
具体的，在本实施例中，所述Q11具体包括：
基于预先获取的MEG系统所采集到的被试者数据，采用公式（2）获取预先获取的MEG系统所采集到的被试者数据X（t）的离散样本X（r）的傅里叶系数。
所述公式（2）为：
Figure 654060DEST_PATH_IMAGE010
；
Figure 42316DEST_PATH_IMAGE011
为MEG系统所采集到的被试者数据X（t）的离散样本X（r）的傅里叶系数。
R表示预先获取的MEG系统所采集到的被试者数据时域采样点数；r表示离散样本X （r）中的第r个样本。
基于所述X（r）的傅里叶系数
Figure 533341DEST_PATH_IMAGE011
，采用公式（3）获取
Figure 919323DEST_PATH_IMAGE012
的具体值，并将所述
Figure 745065DEST_PATH_IMAGE012
的具 体值作为C a 的具体值。
所述公式（3）为：
Figure 937012DEST_PATH_IMAGE013
；
Figure 16963DEST_PATH_IMAGE012
为C a 的估计值。
H表示对矩阵做共轭转置；
Figure 839426DEST_PATH_IMAGE014
为
Figure 637618DEST_PATH_IMAGE011
的共轭转置。
其中，n a 为预先划分MEG系统所采集到的被试者数据的频带中的第a个频带内可用的傅里叶系数数目，并且，
Figure 383988DEST_PATH_IMAGE015
；
其中，#表示计数。
I a 为预先划分的频带中的第a个频带。
k表示总的傅里叶系数数目。
具体的，在本实施例中，所述Q12具体包括：
基于预先获取的MEG系统所采集到空房间中的背景噪声信号，采用公式（4）获取预先获取的MEG系统所采集到的背景噪声信号N（t）的离散样本N（w）的傅里叶系数。
所述公式（4）为：
Figure 521708DEST_PATH_IMAGE016
；
Figure 311810DEST_PATH_IMAGE017
为MEG系统所采集到的背景噪声信号N（t）的离散样本N（w）的傅里叶系数。
W表示MEG系统所采集到的背景噪声信号时域采样点数。
w表示离散样本N（w）中的第w个样本。
基于所述N（w）的傅里叶系数
Figure 862877DEST_PATH_IMAGE017
，采用公式（5）获取
Figure 865468DEST_PATH_IMAGE018
的具体值，并将所述
Figure 160490DEST_PATH_IMAGE018
的具 体值作为B a 的具体值。
所述公式（5）为：
Figure 387072DEST_PATH_IMAGE019
；
Figure 97539DEST_PATH_IMAGE018
为B a 的估计值。
Figure 169400DEST_PATH_IMAGE020
为
Figure 609609DEST_PATH_IMAGE017
的共轭转置。
其中，z a 为预先划分MEG系统所采集到的背景噪声信号的频带中的第a个频带内可用的傅里叶系数数目，并且，
Figure 226666DEST_PATH_IMAGE021
。
具体的，在本实施例中，所述Q13具体包括：
基于KL散度构造频谱匹配标准，所述频谱匹配标准的形式为：
Figure 690008DEST_PATH_IMAGE022
其中，m为频带数目。
Figure 299981DEST_PATH_IMAGE023
表示
Figure 594696DEST_PATH_IMAGE024
和
Figure 631923DEST_PATH_IMAGE025
的KL散度。
基于所述频谱匹配标准、所述C a 参数的具体值和所述B a 参数的具体值，采用期望最大化算法获取单个频带内的源功率A a 参数的具体值和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值。
本实施例中的一种频谱匹配独立成分分析方法，由于通过获取的MEG系统所采集到的被试者数据中单个频带内的噪声功率B a 参数的具体值、单个频带内的源功率A a 参数的具体值、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值，得到该MEG系统所采集到的被试者数据中预先设定的频带内的进行滤波的信号矩阵X a 中的源信号，使得源信号恢复更为准确。
本领域内的技术人员应明白，本发明的实施例可提供为方法、系统或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例，或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质（包括但不限于磁盘存储器、CD-ROM、光学存储器等）上实施的计算机程序产品的形式。
本发明是参照根据本发明实施例的方法、设备（系统）和计算机程序产品的流程图和／或方框图来描述的。应理解可由计算机程序指令实现流程图和／或方框图中的每一流程和／或方框、以及流程图和／或方框图中的流程和／或方框的结合。
应当注意的是，在权利要求中，不应将位于括号之间的任何附图标记理解成对权利要求的限制。词语“包含”不排除存在未列在权利要求中的部件或步骤。位于部件之前的词语“一”或“一个”不排除存在多个这样的部件。本发明可以借助于包括有若干不同部件的硬件以及借助于适当编程的计算机来实现。在列举了若干装置的权利要求中，这些装置中的若干个可以是通过同一个硬件来具体体现。词语第一、第二、第三等的使用，仅是为了表述方便，而不表示任何顺序。可将这些词语理解为部件名称的一部分。
此外，需要说明的是，在本说明书的描述中，术语“一个实施例”、“一些实施例”、“实施例”、“示例”、“具体示例”或“一些示例”等的描述，是指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。
尽管已描述了本发明的优选实施例，但本领域的技术人员在得知了基本创造性概念后，则可对这些实施例作出另外的变更和修改。所以，权利要求应该解释为包括优选实施例以及落入本发明范围的所有变更和修改。
显然，本领域的技术人员可以对本发明进行各种修改和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也应该包含这些修改和变型在内。

Claims (10)

1.一种频谱匹配独立成分分析方法，其特征在于，包括：
S1、获取MEG系统所采集到的被试者数据中单个频带内的噪声功率B a 参数的具体值、单个频带内的源功率A a 参数的具体值、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值；
所述噪声功率B a 参数为MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵在预先划分的频带中的第a个频带内的对应平均值；
所述源功率A a 参数为MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵在预先划分的频带中的第a个频带内的对应平均值；
S2、基于所述MEG系统所采集到的被试者数据中单个频带内的噪声功率B a 参数的具体值、单个频带内的源功率A a 参数的具体值、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值，获取该MEG系统所采集到的被试者数据中预先设定的频带内的进行滤波的信号矩阵X a 中的源信号。
2.根据权利要求1所述的方法，其特征在于，所述S2具体包括：
基于所述MEG系统所采集到的被试者数据中单个频带内的噪声功率B a 参数的具体值、单个频带内的源功率A a 参数的具体值、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值，根据公式（1）计算得到该MEG系统所采集到的被试者数据中第a频带范围上的进行滤波的信号矩阵X a 中的源信号；
所述公式（1）为：
Figure 219786DEST_PATH_IMAGE001
其中，
Figure 930253DEST_PATH_IMAGE002
为该MEG系统所采集到的被试者数据中预先划分的频带中的第a个频带内的进 行滤波的信号矩阵X a 对应恢复的源信号矩阵；
Figure 533273DEST_PATH_IMAGE003
为MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所 采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的转置；
Figure 645585DEST_PATH_IMAGE004
为所述MEG系统所采集到的被试者数据中单个频带内的噪声功率B a 参数的逆；
Figure 246331DEST_PATH_IMAGE005
为单个频带内的源功率A a 参数的逆。
3.根据权利要求2所述的方法，其特征在于，在S1之前还包括：
S0、基于预先获取的MEG系统所采集到的被试者数据、预先获取的背景噪音信号，建立频谱模型，所述频谱模型中包括噪声功率B a 参数、单个频带内的源功率A a 参数、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M参数；
相应的，所述S1具体包括：
基于预先获取的MEG系统所采集到的被试者数据、预先获取的MEG系统采集的空房间中的背景噪声信号和所述频谱模型，获取MEG系统所采集到的被试者数据中单个频带内的噪声功率B a 参数的具体值、单个频带内的源功率A a 参数的具体值、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值。
4.根据权利要求3所述的方法，其特征在于，
所述频谱模型为：
Figure 568728DEST_PATH_IMAGE006
其中，C a 为C（f）在所划分的单个频带内的对应平均值；
C（f）为X（t）的自协方差矩阵的傅里叶变换；
其中，所述X（t）为p维向量，表示从具有p个探头传感器的MEG系统所采集到的被试者数据；
其中，
Figure 850805DEST_PATH_IMAGE007
其中，N（t）为MEG系统所采集到的被试者数据时域内的背景噪声项且为p维向量；
Figure 440793DEST_PATH_IMAGE008
表示MEG系统所采集到的被试者数据时域内的q个独立源且为q维向量；
M是一个
Figure 478019DEST_PATH_IMAGE009
的混合矩阵，表示假设通过该矩阵源信号与噪声信号实现加性混合。
5.根据权利要求4所述的方法，其特征在于，所述S1具体包括：
S11、基于预先获取的MEG系统所采集到的被试者数据，获取C a 参数的具体值；
S12、基于预先获取的背景噪声信号，获取B a 参数的具体值；
S13、基于所述C a 参数的具体值和所述B a 参数的具体值，获取单个频带内的源功率A a 参数的具体值、MEG系统所采集到的被试者数据频域内的噪声项的频谱协方差矩阵和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值。
6.根据权利要求5所述的方法，其特征在于，所述S11具体包括：
基于预先获取的MEG系统所采集到的被试者数据，采用公式（2）获取预先获取的MEG系统所采集到的被试者数据X（t）的离散样本X（r）的傅里叶系数；
所述公式（2）为：
Figure 897499DEST_PATH_IMAGE011
；
Figure 107900DEST_PATH_IMAGE012
为MEG系统所采集到的被试者数据X（t）的离散样本X（r）的傅里叶系数；
R表示预先获取的MEG系统所采集到的被试者数据时域采样点数；r表示离散样本X（r）中的第r个样本；
基于所述X（r）的傅里叶系数
Figure 929226DEST_PATH_IMAGE012
，采用公式（3）获取
Figure 465249DEST_PATH_IMAGE013
的具体值，并将所述
Figure 637605DEST_PATH_IMAGE013
的具体值 作为C a 的具体值；
所述公式（3）为：
Figure 651697DEST_PATH_IMAGE015
；
Figure 593108DEST_PATH_IMAGE013
为C a 的估计值；
H表示对矩阵做共轭转置；
Figure 535919DEST_PATH_IMAGE016
为
Figure 195570DEST_PATH_IMAGE012
的共轭转置；
其中，n a 为预先划分MEG系统所采集到的被试者数据的频带中的第a个频带内可用的傅里叶系数数目，并且，
Figure 482195DEST_PATH_IMAGE017
；
其中，#表示计数；
I a 为预先划分的频带中的第a个频带；
k表示总的傅里叶系数数目。
7.根据权利要求6所述的方法，其特征在于，所述S12具体包括：
基于预先获取的MEG系统所采集到空房间中的背景噪声信号，采用公式（4）获取预先获取的MEG系统所采集到的背景噪声信号N（t）的离散样本N（w）的傅里叶系数；
所述公式（4）为：
Figure 137167DEST_PATH_IMAGE019
；
Figure 890360DEST_PATH_IMAGE020
为MEG系统所采集到的背景噪声信号N（t）的离散样本N（w）的傅里叶系数；
W表示MEG系统所采集到的背景噪声信号时域采样点数；
w表示离散样本N（w）中的第w个样本；
基于所述N（w）的傅里叶系数
Figure 896362DEST_PATH_IMAGE020
，采用公式（5）获取
Figure 393202DEST_PATH_IMAGE021
的具体值，并将所述
Figure 481111DEST_PATH_IMAGE021
的具体值 作为B a 的具体值；
所述公式（5）为：
Figure 139625DEST_PATH_IMAGE023
；
Figure 305028DEST_PATH_IMAGE021
为B a 的估计值；
Figure DEST_PATH_IMAGE024
为
Figure 667876DEST_PATH_IMAGE020
的共轭转置；
其中，z a 为预先划分MEG系统所采集到的背景噪声信号的频带中的第a个频带内可用的傅里叶系数数目，并且，
Figure DEST_PATH_IMAGE026
。
8.根据权利要求7所述的方法，其特征在于，所述S13具体包括：
基于KL散度构造频谱匹配标准，所述频谱匹配标准的形式为：
Figure DEST_PATH_IMAGE028
其中，m为频带数目；
Figure DEST_PATH_IMAGE030
表示
Figure DEST_PATH_IMAGE032
和
Figure DEST_PATH_IMAGE034
的KL散度；
基于所述频谱匹配标准、所述C a 参数的具体值和所述B a 参数的具体值，采用期望最大化算法获取单个频带内的源功率A a 参数的具体值和MEG系统所采集到的被试者数据频域内q个独立源频谱协方差矩阵的加性混合矩阵M的具体值。
9.根据权利要求8所述的方法，其特征在于，
所述MEG系统所采集到的被试者数据为用户处于背景剩磁消减后的磁屏蔽房内并在相同声音信号的刺激下静止坐下、静止站立以及身体和头部移动时所采集的MEG被试者数据；
所述背景剩磁消减后的磁屏蔽房内的背景磁场低于50nT；
预先获取的背景噪音信号为背景剩磁消减后的磁屏蔽房内的预先设定时间段内的空的磁屏蔽房内的磁信号。
10. 一种频谱匹配独立成分分析系统，其特征在于，所述系统包括：
至少一个处理器；以及
与所述处理器通信连接的至少一个存储器，其中，所述存储器存储有可被所述处理器执行的程序指令，所述处理器调用所述程序指令能够执行如权利要求1至9中任一所述的频谱匹配独立成分分析方法。
CN202111583666.0A 2021-12-22 2021-12-22 一种频谱匹配独立成分分析方法及系统 Active CN113951887B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111583666.0A CN113951887B (zh) 2021-12-22 2021-12-22 一种频谱匹配独立成分分析方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111583666.0A CN113951887B (zh) 2021-12-22 2021-12-22 一种频谱匹配独立成分分析方法及系统

Publications (2)

Publication Number Publication Date
CN113951887A true CN113951887A (zh) 2022-01-21
CN113951887B CN113951887B (zh) 2022-03-25

Family

ID=79473678

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111583666.0A Active CN113951887B (zh) 2021-12-22 2021-12-22 一种频谱匹配独立成分分析方法及系统

Country Status (1)

Country Link
CN (1) CN113951887B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115018018A (zh) * 2022-08-05 2022-09-06 北京航空航天大学杭州创新研究院 一种抑制背景噪声的双重空间滤波方法

Citations (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6544170B1 (en) * 1999-06-21 2003-04-08 Shimadzu Corporation Biosignal measuring method and apparatus
US20050234329A1 (en) * 2004-04-15 2005-10-20 Kraus Robert H Jr Noise cancellation in magnetoencephalography and electroencephalography with isolated reference sensors
WO2006072150A1 (en) * 2005-01-07 2006-07-13 K.U. Leuven Research And Development Muscle artifact removal from encephalograms
US20140107464A1 (en) * 2012-06-29 2014-04-17 Commissariat A L'energie Atomique Et Aux Energies Alternatives Method for Locating a Brain Activity
CN105147288A (zh) * 2015-07-23 2015-12-16 中国科学院苏州生物医学工程技术研究所 脑磁源强度定位方法
US20170311832A1 (en) * 2014-11-13 2017-11-02 Mensia Technologies Scoring method based on improved signals analysis
CN107562188A (zh) * 2017-07-17 2018-01-09 北京大学 基于原子磁强计的脑机接口系统及其使用方法
CN109243483A (zh) * 2018-10-17 2019-01-18 西安交通大学 一种含噪频域卷积盲源分离方法
CN111543949A (zh) * 2020-05-13 2020-08-18 北京航空航天大学 一种基于脑磁图与脑电图的儿童asd诊断装置
CN111667846A (zh) * 2020-05-14 2020-09-15 北京声智科技有限公司 一种盲源分离方法及装置
CN112674773A (zh) * 2020-12-22 2021-04-20 北京航空航天大学 基于Tucker分解和ripple时间窗的脑磁图源定位方法和装置
CN112914578A (zh) * 2021-01-20 2021-06-08 季华实验室 Meg源定位方法及系统
US20210186400A1 (en) * 2018-08-24 2021-06-24 University Of Washington Apparatuses, systems, and methods for suppression of artifacts in non-invasive electromagnetic recordings

Patent Citations (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6544170B1 (en) * 1999-06-21 2003-04-08 Shimadzu Corporation Biosignal measuring method and apparatus
US20050234329A1 (en) * 2004-04-15 2005-10-20 Kraus Robert H Jr Noise cancellation in magnetoencephalography and electroencephalography with isolated reference sensors
WO2006072150A1 (en) * 2005-01-07 2006-07-13 K.U. Leuven Research And Development Muscle artifact removal from encephalograms
US20140107464A1 (en) * 2012-06-29 2014-04-17 Commissariat A L'energie Atomique Et Aux Energies Alternatives Method for Locating a Brain Activity
US20170311832A1 (en) * 2014-11-13 2017-11-02 Mensia Technologies Scoring method based on improved signals analysis
CN105147288A (zh) * 2015-07-23 2015-12-16 中国科学院苏州生物医学工程技术研究所 脑磁源强度定位方法
CN107562188A (zh) * 2017-07-17 2018-01-09 北京大学 基于原子磁强计的脑机接口系统及其使用方法
US20210186400A1 (en) * 2018-08-24 2021-06-24 University Of Washington Apparatuses, systems, and methods for suppression of artifacts in non-invasive electromagnetic recordings
CN109243483A (zh) * 2018-10-17 2019-01-18 西安交通大学 一种含噪频域卷积盲源分离方法
CN111543949A (zh) * 2020-05-13 2020-08-18 北京航空航天大学 一种基于脑磁图与脑电图的儿童asd诊断装置
CN111667846A (zh) * 2020-05-14 2020-09-15 北京声智科技有限公司 一种盲源分离方法及装置
CN112674773A (zh) * 2020-12-22 2021-04-20 北京航空航天大学 基于Tucker分解和ripple时间窗的脑磁图源定位方法和装置
CN112914578A (zh) * 2021-01-20 2021-06-08 季华实验室 Meg源定位方法及系统

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
ESCUDERO J ET AL.: "Artifact removal in magnetoencephalogram background activity with independent component analysis", 《IEEE TRANS BIOMED ENG》 *
武振华等: "独立成分分析在生物医学信号处理中的应用", 《国外医学(生物医学工程分册)》 *
王斌等: "基于独立元分析的MEG（脑磁图）数据分析和处理", 《生物物理学报》 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115018018A (zh) * 2022-08-05 2022-09-06 北京航空航天大学杭州创新研究院 一种抑制背景噪声的双重空间滤波方法

Also Published As

Publication number Publication date
CN113951887B (zh) 2022-03-25

Similar Documents

Publication Publication Date Title
Taulu et al. Applications of the signal space separation method
Taulu et al. Presentation of electromagnetic multichannel data: the signal space separation method
Unser et al. A review of wavelets in biomedical applications
US11083401B2 (en) Electric field encephalography: electric field based brain signal detection and monitoring
Moiseev et al. Application of multi-source minimum variance beamformers for reconstruction of correlated neural activity
Kobayashi et al. Principal component elimination method for the improvement of S/N in evoked neuromagnetic field measurements
CN110353673B (zh) 一种基于标准互信息的脑电通道选择方法
CN106618497A (zh) 复杂环境下基于信道状态信息监测睡眠的方法
WO1994028789A1 (en) Apparatus and process for determining the sources of biomagnetic activity
Vigário et al. Independence: A new criterion for the analysis of the electromagnetic fields in the global brain?
CN112914578B (zh) Meg源定位方法及系统
TW201416694A (zh) 一種重建多通道核磁共振影像的方法
CN113951887B (zh) 一种频谱匹配独立成分分析方法及系统
Boux et al. Dictionary-free MR fingerprinting parameter estimation via inverse regression
Sekihara et al. MEG spatio-temporal analysis using a covariance matrix calculated from nonaveraged multiple-epoch data
Taulu et al. The signal space separation method
Sekihara et al. Estimating neural sources from each time-frequency component of magnetoencephalographic data
Arvinti et al. Performance of Fourier versus Wavelet analysis for magnetocardiograms using a SQUID-acquisition system
Kawakatsu Application of ICA to MEG noise reduction
WO2000010454A1 (en) Functional brain imaging from magnetoencephalographic data
CN113171098A (zh) 用于大脑晃动中脑磁信号的检测装置和检测方法
CN115018018B (zh) 一种抑制背景噪声的双重空间滤波方法
Taulu Processing of weak magnetic multichannel signals: the signal space separation method
Mino et al. Neuromagnetic source analysis with a 64-channel SQUID system and MR imaging
Arvinti-Costache et al. A wavelet based baseline drift correction method for fetal magnetocardiograms

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant