CN103190898A - 心磁信号噪声自适应滤波消除设计方法 - Google Patents

心磁信号噪声自适应滤波消除设计方法 Download PDF

Info

Publication number
CN103190898A
CN103190898A CN2013101425659A CN201310142565A CN103190898A CN 103190898 A CN103190898 A CN 103190898A CN 2013101425659 A CN2013101425659 A CN 2013101425659A CN 201310142565 A CN201310142565 A CN 201310142565A CN 103190898 A CN103190898 A CN 103190898A
Authority
CN
China
Prior art keywords
time
frequency
singular value
conversion
matrix
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
CN2013101425659A
Other languages
English (en)
Other versions
CN103190898B (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.)
Ningbo Lidou Intelligent Technology Co ltd
Original Assignee
Individual
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 Individual filed Critical Individual
Priority to CN201310142565.9A priority Critical patent/CN103190898B/zh
Publication of CN103190898A publication Critical patent/CN103190898A/zh
Application granted granted Critical
Publication of CN103190898B publication Critical patent/CN103190898B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)

Abstract

心磁信号噪声自适应滤波消除设计方法，包括以下步骤：(1)对实时采集到的心磁信号进行S变换，获得信号时频域内特征矩阵S变换时频矩阵；(2)在S变换中引入时频调节因子
Figure 2013101425659100004DEST_PATH_IMAGE002
、
Figure DEST_PATH_IMAGE004
构建广义S变换，利用广义S变换时频调节因子调节信号时频特征的时频分辨率；(3)采用奇异值分解方法分解广义S变换二维时频特征矩阵；(4)计算单个奇异值时域信号与噪声的相关系数，取相关系数r<0.5且单奇异值占总奇异值比例大于10%奇异值所对应的矩阵区域为有效心磁信号奇异值分布区域，对有效心磁信号时频域进行反S变换重构信号。本发明与现有技术相比，采用广义S变换和奇异值分解方法，可以在无参考噪声数据情形下解决背景噪声滤除问题。

Description

心磁信号噪声自适应滤波消除设计方法
技术领域
本发明属于电磁信号处理领域，涉及一种基于广义S变换和奇异值分解的心磁信号噪声自适应滤波消除设计方法。
背景技术
心脏疾病是危害人类健康的主要疾病之一。目前，心脏功能检测和诊断方法主要采用心电图来检测和分析心脏电活动。随着磁传感器件的发展，高温射频超导量子干涉仪作为灵敏度极高的磁传感器件广泛的应用于心脏信号的测量中。高灵敏度的高温射频超导量子干涉仪在无屏蔽室内采集到的心磁信号中存在大量的背景噪声，其主要成份为50HZ的市电干扰及其它环境噪声。如何从包含大量背景噪声的含噪信号中有效的提取心磁信号及相关特征信息是心脏病诊断与治疗的关键。在现已报道的心磁信号噪声自适应滤波消除设计方法中，主要包括三种方法：一种是基于自适应理论的背景噪声消除方法，自适应理论算法适用于心磁信号及环境噪声可同步采集情形，但不能解决无参考噪声数据情形下的背景噪声滤除问题。另一种是基于SVD分解和自适应滤波的心磁信号消噪方法。该方法对无参考噪声情形下的近周期心磁信号消噪问题具有较好消除效果，但主要缺点是其所分解的Hankel矩阵图像特征不能提取信号及背景噪声基本频域特征，且在消除50HZ市电干扰时，仍需采用自适应滤波方法。第三种是利用小波变换来处理心磁信号，采用多种小波基函数在不同尺度下进行消噪研究。结果显示，采用symlet8小波函数进行消噪能取得了较好的滤波效果，其局局限性在于选取合适的小波基函数和尺度因子与背景噪声强度有关，选取过程较为复杂，且在高频段时需要强制消噪处理。
发明内容
本发明要解决的技术问题是，克服现有心磁信号滤波设计方法中存在的上述不足和满足非平稳信号处理未来发展的需要，提供了一种基于广义S变换和奇异值分解的心磁信号噪声自适应滤波消除设计方法。该方法通用性强，可以在无参考噪声数据情形下解决背景噪声滤除问题；在消除50HZ市电干扰及其它定频干扰时，无需采用其它自适应滤波方法；该方法实现简单，噪声抑制比高，运行速度快，能在较高信噪比的条件下，利用较少奇异值就能取得良好的滤波效果。
本发明解决其技术问题所采用的技术方案是：
心磁信号噪声自适应滤波消除设计方法，包括以下步骤： (1) 对实时采集到的心磁信号进行S变换，获得信号时频域内特征矩阵S变换时频矩阵； (2) 根据S变换时频矩阵时频分辨率的实时要求高低，在S变换中引入时频调节因子                                                
Figure 889279DEST_PATH_IMAGE001
、
Figure 2013101425659100002DEST_PATH_IMAGE002
构建广义S变换，利用广义S变换时频调节因子调节信号时频特征的时频分辨率； (3) 采用奇异值分解方法分解广义S变换二维时频特征矩阵，得到左特征时域矩阵、右特征频域矩阵及特征值组成的对角阵，将单个奇异值对角阵分别与左右特征时域矩阵相乘得到单奇异值S变换时频矩阵，单奇异值S变换时频矩阵经反S变换即可得到单个主要奇异值所对应的信号时域成份；(4) 计算单个奇异值时域信号与高斯白噪声（均值为0，方差为0.1）的相关系数，取相关系数r<0.5且单个奇异值占总奇异值比例大于10%的奇异值所对应奇异矩阵区域为有效心磁信号奇异值分布区域，对有效心磁信号时频域进行反S变换重构信号从而实现心磁信号自适应时频滤波。
所述步骤（1）中，所述心磁信号是指工作在无磁屏蔽室条件下利用高温射频超导量子干涉仪系统所测量到的心脏电心理过程在躯干上面产生的时变磁场强度信号。
所述步骤（1）中，所述S变换是指以高斯函数为窗函数的信号短时傅里叶变换，其表达式为
Figure 619468DEST_PATH_IMAGE003
，式中，
Figure 2013101425659100002DEST_PATH_IMAGE004
为时间，
Figure 274572DEST_PATH_IMAGE005
为频率，
Figure 2013101425659100002DEST_PATH_IMAGE006
为时移因子，
Figure 607464DEST_PATH_IMAGE007
为虚数，
Figure 2013101425659100002DEST_PATH_IMAGE008
为表示实时连续心磁信号，
Figure 317710DEST_PATH_IMAGE009
为高斯窗函数，其中。
所述步骤（2）中，所述广义S变换是指在S变换中引入时频调节因子
Figure 77856DEST_PATH_IMAGE011
、
Figure 282572DEST_PATH_IMAGE002
，其表达式为
Figure 2013101425659100002DEST_PATH_IMAGE012
。其中，
Figure 419156DEST_PATH_IMAGE001
为高斯窗幅度拉伸因子，
Figure 25717DEST_PATH_IMAGE002
为频率尺度拉伸因子，
Figure 956764DEST_PATH_IMAGE004
为时间，
Figure 711094DEST_PATH_IMAGE005
为频率，
Figure 385789DEST_PATH_IMAGE006
为时移因子，
Figure 846857DEST_PATH_IMAGE007
为虚数，
Figure 211455DEST_PATH_IMAGE008
为表示实时连续心磁信号。当且仅当
Figure 390763DEST_PATH_IMAGE013
，且
Figure 2013101425659100002DEST_PATH_IMAGE014
时，广义S变换为StockWell所提出的
Figure 603570DEST_PATH_IMAGE015
变换。当
Figure 2013101425659100002DEST_PATH_IMAGE016
或时，频率分辨率提高，时间分辨率下降；当
Figure 2013101425659100002DEST_PATH_IMAGE018
或时，时间分辨率提高，频率分辨率下降。
所述步骤（3）中，所述奇异值分解方法是指一种非退化正交矩阵分解法，在本发明中是指心磁信号经广义S变换后的二维时频特征矩阵正交分解方法。
本发明的原理是：首先，对实时采集到的心磁信号进行S变换，获得时频域信号特征。然后，利用广义S变换调节时频特征的时频分辨率，使信号特征与噪声特征尽可能的分离。最后，采用奇异值分解方法分解时频特征，提取单个主要奇异值信号时域成份，利用单个奇异值时域信号与噪声的相关系数确定有效心磁信号奇异值分布区域，通过保留奇异值较大且相关度较小的奇异值重构信号，从而实现心磁信号自适应时频滤波，达到滤除噪声提取心磁信号的目的。
1、S变换及广义S变换原理
S变换起源于短时傅里叶变换，也可由小波变换引出。设一维连续信号的小波变换定义如下
Figure 2013101425659100002DEST_PATH_IMAGE020
                   (1) ，
式中，
Figure 13112DEST_PATH_IMAGE006
为时移因子，
Figure 448773DEST_PATH_IMAGE021
为尺度因子，为小波基函数，
Figure 158103DEST_PATH_IMAGE023
为信号
Figure 46424DEST_PATH_IMAGE008
在尺度
Figure 804296DEST_PATH_IMAGE021
下的小波变换系数。当小波变换基函数选用高斯窗函数
Figure 2013101425659100002DEST_PATH_IMAGE024
时，其对应的小波变换为
Figure 32146DEST_PATH_IMAGE025
                    (2)，
式中，
Figure 2013101425659100002DEST_PATH_IMAGE026
为信号
Figure 722497DEST_PATH_IMAGE008
在高斯窗基函数下的小波变换系数，尺度
Figure 973481DEST_PATH_IMAGE027
。对(2)式两边同时乘以相位修正因子
Figure 2013101425659100002DEST_PATH_IMAGE028
，取，并对其幅值进行
Figure 2013101425659100002DEST_PATH_IMAGE030
修正可得
Figure 676787DEST_PATH_IMAGE015
变换为
Figure 727920DEST_PATH_IMAGE031
                    (3)，
式中，
Figure 590834DEST_PATH_IMAGE006
为时移因子，
Figure 815142DEST_PATH_IMAGE005
为频率，
Figure 2013101425659100002DEST_PATH_IMAGE032
为信号
Figure 562124DEST_PATH_IMAGE008
在不同时间不同频率下经S变换后的二维时频特征矩阵。
为了有效利用快速FT算法，将(3)式重写可得变换快速算法为
Figure 134368DEST_PATH_IMAGE033
                       (4)，
式中，
Figure 2013101425659100002DEST_PATH_IMAGE034
为频移因子，
Figure 37733DEST_PATH_IMAGE035
为信号
Figure 156999DEST_PATH_IMAGE008
的傅里叶变换，
Figure 2013101425659100002DEST_PATH_IMAGE036
为信号
Figure 284355DEST_PATH_IMAGE008
经快速S变换后得到的二维时频特征矩阵。
在离散情况下，信号
Figure 121861DEST_PATH_IMAGE008
经采样后得到
Figure 622725DEST_PATH_IMAGE037
点离散序列
Figure 2013101425659100002DEST_PATH_IMAGE038
，
Figure 799760DEST_PATH_IMAGE039
为采样周期， 则直接计算离散的
Figure 35700DEST_PATH_IMAGE015
变换
Figure 422819DEST_PATH_IMAGE041
为
Figure 2013101425659100002DEST_PATH_IMAGE042
                 (5)，
式中，
Figure 667987DEST_PATH_IMAGE043
为离散采样序列，
Figure 558582DEST_PATH_IMAGE037
为采样点数，
Figure 2013101425659100002DEST_PATH_IMAGE044
为时间离散采样序列，
Figure 228073DEST_PATH_IMAGE045
为频率离散采样序列。
Figure 836909DEST_PATH_IMAGE041
为离散普通S变换得到的二维时频特征矩阵。
采用快速FT算法计算
Figure 948085DEST_PATH_IMAGE015
变换
Figure 2013101425659100002DEST_PATH_IMAGE046
为
Figure 896449DEST_PATH_IMAGE047
                       (6)，
式中，
Figure 2013101425659100002DEST_PATH_IMAGE048
为离散频移因子，
Figure 802088DEST_PATH_IMAGE049
为信号序列
Figure 835903DEST_PATH_IMAGE043
的离散傅里叶变换。其中，
Figure 2013101425659100002DEST_PATH_IMAGE050
，
Figure 485191DEST_PATH_IMAGE051
为
Figure 350378DEST_PATH_IMAGE049
在频域上的扩维。
由于
Figure 426919DEST_PATH_IMAGE015
变换的高斯窗函数满足
Figure 2013101425659100002DEST_PATH_IMAGE052
                             （7），
式中，高斯窗函数为
Figure 948030DEST_PATH_IMAGE053
。因此，
Figure 155937DEST_PATH_IMAGE015
变换后的函数对时间的连续积分可以写成
   
Figure 2013101425659100002DEST_PATH_IMAGE054
              (8)，
式中，
Figure 78893DEST_PATH_IMAGE055
为信号
Figure 388652DEST_PATH_IMAGE008
的连续S变换，
Figure 397059DEST_PATH_IMAGE035
为信号
Figure 450466DEST_PATH_IMAGE008
的连续傅里叶变换。由（8）式可知，反变换可由
Figure 646272DEST_PATH_IMAGE035
反FT即可得
Figure 938713DEST_PATH_IMAGE008
，即
                         (9)，
当噪声中存在频率成份与信号成份相近的情况时，采用
Figure 733494DEST_PATH_IMAGE015
变换二维等高图波峰去辨识信号频率成份时，相近时间或频率波峰主成份重叠在一起，从而影响信号的参数提取。为了有效的解决变换的时频分辨率的可调性，有效提取信号的主时频区域，在S变换中引入时频调节因子
Figure 17025DEST_PATH_IMAGE001
，
Figure 24DEST_PATH_IMAGE002
。因此，广义的
Figure 395233DEST_PATH_IMAGE015
变换定义为
Figure 878780DEST_PATH_IMAGE012
                   (10)，
式中，
Figure 904505DEST_PATH_IMAGE001
为高斯窗幅度拉伸因子，
Figure 171538DEST_PATH_IMAGE002
为频率尺度拉伸因子。当且仅当
Figure 308121DEST_PATH_IMAGE013
，且
Figure 711421DEST_PATH_IMAGE014
时，广义S变换为StockWell所提出的
Figure 908047DEST_PATH_IMAGE015
变换；当
Figure 662376DEST_PATH_IMAGE016
或
Figure 337071DEST_PATH_IMAGE017
时，频率分辨率提高时间分辨率下降，当
Figure 798140DEST_PATH_IMAGE018
或
Figure 227984DEST_PATH_IMAGE019
时，时间分辨率提高频率分辨率下降。由于
Figure 141713DEST_PATH_IMAGE001
，
Figure 682416DEST_PATH_IMAGE002
均为实数，反变换对窗函数进行积分时仍然满足(8)式，因此，广义
Figure 263570DEST_PATH_IMAGE015
变换反变换与反变换表达式相同，均为(9)式。
由于S变换后的二维时频特征矩阵表示对应该时间点的频率信息，矩阵模值的大小反映了信号在该时频域内所含分量的大小，因此，滤波的过程就是要尽可能多保留信号在时频域的分量，在信号所对应的分量前乘以较大权值，而在噪声所对应的时频分量前乘以较小权值。S变换时频滤波模型可表示为
                          (11)，
式中，
Figure DEST_PATH_IMAGE058
为自适应滤波器，
Figure 278909DEST_PATH_IMAGE059
为滤波后的时频矩阵。在阀值滤波情况下，信号时频域通过区域为
Figure DEST_PATH_IMAGE060
，阻带通过区域为
Figure 714570DEST_PATH_IMAGE061
。
2、时频域特征奇异值分解原理
广义S变换后的时频矩阵为
Figure DEST_PATH_IMAGE062
的方阵，由矩阵的奇异值与特征向量的关系可得
Figure 158320DEST_PATH_IMAGE063
                                (12)，
式中，
Figure 312221DEST_PATH_IMAGE015
为离散广义S变换二维时频特征矩阵，
Figure DEST_PATH_IMAGE064
为广义S变换时频矩阵的第个奇异值，
Figure 422577DEST_PATH_IMAGE065
，
Figure DEST_PATH_IMAGE066
，为相应的左右特征向量，分别对应广义S变换时频矩阵的时域特征和频域特征。由式(12)可知，信号奇异值的大小对应着信号各时域或频域主要成份在原信号中所占比重的大小，奇异值越大，其对应的时域向量越接近原信号的主要成分，奇异值较小且能量分布较分散的越趋向于噪声。广义S变换后的二维时频特征矩阵奇异值分解可描述为
Figure DEST_PATH_IMAGE068
                           （13），
式中，
Figure 675496DEST_PATH_IMAGE069
，
Figure DEST_PATH_IMAGE070
分别为左特征向量
Figure 33796DEST_PATH_IMAGE066
组成的矩阵和右特征向量所组成的矩阵，
Figure 291919DEST_PATH_IMAGE071
为矩阵
Figure 154833DEST_PATH_IMAGE070
的共轭转置。左特征矩阵由一组代表信号时域信息的酉向量构成，右特征矩阵由一组代表频域信息的酉向量构成。
Figure DEST_PATH_IMAGE072
为一组由奇异值组成的对角阵，其中特征值。由于奇异值的大小反应对应信号主要成份的大小，较大的奇异值主要为信号的主成份，较小的奇异值对应信号的噪声分量，因此，滤波的过程就是保留对角阵中较大的奇异值舍去较小的奇异值重构实现滤波。即
Figure DEST_PATH_IMAGE074
，其中
Figure 581583DEST_PATH_IMAGE075
为二维时频特征矩阵
Figure 803617DEST_PATH_IMAGE015
的轶。
Figure 216144DEST_PATH_IMAGE015
变换具有与小波变换类似的多分辩特性，利用S变换多分辩时频特性，可以构建反S变换时频滤波方法。S变换在傅里叶变换的基础上引入与分析信号频率相关的高斯窗函数，在继承傅里叶变换频域分析和小波变换多分辨特性等优点的同时，采用了窗函数宽度随频率呈反比变化的高斯窗函数，在低频段时窗较宽，从而获得较高的频率分辩率；在高频段时时窗较窄，故可获得较高的时间分辨率。因而广泛的应用于非平稳信号分析与处理中。
本发明与现有技术相比，采用基于广义S变换和奇异值分解的心磁信号噪声自适应滤波消除设计方法，可以根据信号的背景噪声强度自适应的提取心磁信号，可以在无参考噪声的情况下实现较好的滤波效果；在消除50HZ市电干扰及其它定频干扰时，无需采用其它自适应滤波方法；该方法实现简单，噪声抑制比高，运行速度快，利用较少奇异值就能取得良好的滤波效果。该设计方法通用性强，适合于无屏蔽无遮掩情况下受高信噪比干扰的心磁信号数据的滤波处理和参数提取。
附图说明
图1为以理想心磁信号受幅值为0.1、频率为50HZ正弦市电干扰和均值为0、方差为0.05的高斯白噪声叠加干扰进行滤波设计的方法流程图； 
图2（a）为受干扰的心磁信号图；
图2（b）为心磁信号S变换三维时频图；
图2（c）为心磁信号S变换二维时频图；
图2（d）为心磁信号广义S变换二维时频图；
图2（e）为滤除50HZ市电干扰后的心磁图；
图2（f）为滤除50HZ市电干扰后的S变换三维时频图；
图3（a）为r=2时频等高图；
图3（b）为r=4时频等高图；
图3（c）为r=8时频等高图；
图3（d）为r=16时频等高图；
图3（e）为r=2时反变换滤波结果；
图3（f）为r=4时反变换滤波结果；
图3（g）为r=8时反变换滤波结果；
图3（h）为r=16时反变换滤波结果；
图4为单个奇异值占总奇异值比例
Figure DEST_PATH_IMAGE076
及占总能量比例变化曲线图；
图5（a）为第一个奇异值对应该的信号域图；
图5（b）为第二个奇异值
Figure 907949DEST_PATH_IMAGE079
对应该的信号域图；
图5（c）为第三个奇异值
Figure DEST_PATH_IMAGE080
对应该的信号域图；
图5（d）为第四个奇异值
Figure 300884DEST_PATH_IMAGE081
对应的信号域图；
图6（a）为滤除50HZ市电干扰后的心磁图曲线；
图6（b）为r=3时广义S变换SVD滤波结果；
图6（c）为r=2时广义S变换SVD滤波结果；
图6（d）为r=1时S变换SVD滤波结果；
图6（e）为去市电干扰后的信号时频图；
图6（f）为r=3时信号时频图；
图6（g）为r=2时信号时频图；
图6（h）为r=1时信号时频图。
具体实施方式
以下结合附图和实施例对本发明作进一步说明。
参照图1，心磁信号噪声自适应滤波消除设计方法，包括以下步骤： (1) 对实时采集到的心磁信号进行S变换，获得信号时频域内特征矩阵S变换时频矩阵；(2) 根据S变换时频矩阵时频分辨率的实时要求高低，在S变换中引入时频调节因子
Figure 138390DEST_PATH_IMAGE001
、
Figure 704501DEST_PATH_IMAGE002
构建广义S变换，利用广义S变换时频调节因子调节信号时频特征的时频分辨率； (3) 采用奇异值分解方法分解广义S变换二维时频特征矩阵，得到左特征时域矩阵、右特征频域矩阵及特征值组成的对角阵，将单个奇异值对角阵分别与左右特征时域矩阵相乘得到单奇异值S变换时频矩阵，单奇异值S变换时频矩阵经反S变换即可得到单个主要奇异值所对应的信号时域成份；(4) 计算单个奇异值时域信号与噪声的相关系数，取相关系数较小且奇异值较大单奇异矩阵区域为有效心磁信号奇异值分布区域，对有效心磁信号时频域进行反S变换重构信号从而实现心磁信号自适应时频滤波。
所述步骤（1）中，所述心磁信号是指工作在无磁屏蔽室条件下利用高温射频超导量子干涉仪系统所测量到的心脏电心理过程在躯干上面产生的时变磁场强度信号。
所述步骤（1）中，所述S变换是指以高斯函数为窗函数的信号短时傅里叶变换，其表达式为
Figure 678273DEST_PATH_IMAGE003
，式中，
Figure 242109DEST_PATH_IMAGE004
为时间，
Figure 566912DEST_PATH_IMAGE005
为频率，为时移因子，
Figure 762062DEST_PATH_IMAGE007
为虚数，
Figure 559117DEST_PATH_IMAGE008
为表示实时连续心磁信号，
Figure 43319DEST_PATH_IMAGE009
为高斯窗函数，其中
Figure 154494DEST_PATH_IMAGE010
。
所述步骤（2）中，所述广义S变换是指在S变换中引入时频调节因子
Figure 165176DEST_PATH_IMAGE011
、
Figure 70815DEST_PATH_IMAGE002
，其表达式为
Figure 104630DEST_PATH_IMAGE012
。其中，
Figure 550655DEST_PATH_IMAGE001
为高斯窗幅度拉伸因子，
Figure 619105DEST_PATH_IMAGE002
为频率尺度拉伸因子，
Figure 757962DEST_PATH_IMAGE004
为时间，
Figure 279073DEST_PATH_IMAGE005
为频率，
Figure 463542DEST_PATH_IMAGE006
为时移因子，为虚数，
Figure 696258DEST_PATH_IMAGE008
为表示实时连续心磁信号。当且仅当
Figure 580031DEST_PATH_IMAGE013
，且
Figure 633438DEST_PATH_IMAGE014
时，广义S变换为StockWell所提出的
Figure 410901DEST_PATH_IMAGE015
变换。当
Figure 829244DEST_PATH_IMAGE016
或
Figure 121685DEST_PATH_IMAGE017
时，频率分辨率提高，时间分辨率下降；当
Figure 916466DEST_PATH_IMAGE018
或
Figure 610752DEST_PATH_IMAGE019
时，时间分辨率提高，频率分辨率下降。
所述步骤（3）中，所述奇异值分解方法是指一种非退化正交矩阵分解法，在本发明中是指心磁信号经广义S变换后的二维时频特征矩阵正交分解方法。
应用实施例：
本发明方法适合任意无磁屏蔽室条件下利用高温射频超导量子干涉仪系统所测量到的心磁信号。为了说明的方便，下面以理想心磁信号受幅值为0.1、频率为50HZ正弦市电干扰和均值为0、方差为0.05的高斯白噪声叠加干扰为例进行滤波设计说明。
设计时，主要包括以下几点：首先，对仿真信号进行采样，得到离散心磁信号，利用快速S变换至二维时频域。根据快速S变换二维时频域内50HZ市电干扰时频分布，调节时频调节因子
Figure 197067DEST_PATH_IMAGE001
，，本设计中
Figure DEST_PATH_IMAGE082
，
Figure 575276DEST_PATH_IMAGE083
。采用广义S变换阀值滤波方法对50HZ市电干扰进行滤波。其次，对滤除市电干扰后的时频域进行奇异值分解，判断各奇异值对应的时域信号是否与心磁信号相关，若相关，则保留相近的奇异值；若不相关，则将奇异值置零。最后通后反广义S变换进行重构信号，实现自适应滤波。
心磁信号经S变换和广义S变换后，信号的时频特征能有效的表征在二维时频图上。图2(a)为仿真心磁信号受幅值为0.1、频率为50HZ正弦市电干扰和均值为0、方差为0.05的高斯噪声干扰信号图，经S变换后的三维时频图如图2(b)所示，从图中可以看出50HZ市电信号成份及对应的时间，心磁信号的主能量相对较大，分布在频率相对较低的频域内，其二维映射时频图如图2(c)所示，从图2(c)可以清晰的看出50HZ频率成份及高频段的噪声分布情况。图2(d)为经广义S变换后的二维时频图，从图2(d)中可以看出，50HZ干扰信号与噪声及心磁信号能相对有效区分，提高了分析信号的时频分辨率。图2(e)为采用广义S变换阀值滤波滤除50HZ市电干扰的心磁信号图，其对应该的三维时频图如图2(f)所示。由图2(f)可知，50HZ市电干扰得到约99%滤除。
滤除50HZ市电干扰后的信号经广义S变换(
Figure DEST_PATH_IMAGE084
，
Figure 999435DEST_PATH_IMAGE085
)得到二维时频特征矩阵，采用奇异值分解方法对其时频矩阵行分解，当奇异值分别取前2个、4个、8个和16个奇异值时，其二维时频等高图分别如图3(a)、图3(b)、图3(c)和图3(d)所示，从二维等高图的变化趋势可知，随着奇异值个数的增加，时频图中噪声分量逐渐增加。经广义反S变换后，其对应的时域信号分别如图3(e)、图3(f)、图3(g)和图3(h)所示，从图3(e)和图3(f)中可以看出，当奇异值个数分别为2、4时具有较好的局部滤波效果，随着所取奇异值个数的增加，还原后的信号中噪声分量逐渐增加。
图4为单个奇异值比例
Figure 87477DEST_PATH_IMAGE076
与单个奇异值能量比例
Figure 292193DEST_PATH_IMAGE077
曲线，由图4可知，第一个奇异值与第二个奇异值占的总奇异值比例和能量比例都相对较大，从第三个奇异值开始，无论是奇异值点比例还是能量比例都相对较小，且比例分布相对变化不大。将前四个奇异分值所对应的时域成份进行提取，其主要特征如图5（a）、图5（b）、图5（c）、图5（d）所示，可知，第一个奇异值和第二个奇异值所对应的信号特征与有遮掩屏蔽理想环境下心磁图P、Q、R、S、T波形特征一致，而第三个奇异值和第四个奇异值所对应的特征分量与噪声分布相似。综合分析信号奇异值分布比例、能量分布比例情况及对应的前四个奇异值的时域特征，可以确定第一个和第二个奇异值对应的是心磁信号特征，而第三个及以后的奇异值主要为噪声特征。因此，滤波时可以保留前两个主奇异值即可实现时频分解滤波。
图6(a)与图6(e)为受背景噪声干扰信号的时域图和广义S变换时频图。由图6(e)可知，信号主能量主要集中在
Figure DEST_PATH_IMAGE086
和附近局部时频域内且频率相对较低，背景噪声主要分布在频率相对较的频带内。对图6(e)所示的时频域进行SVD分解，分别取奇异值个数分别为r=3、r=2和r=1时进行滤波，其对应该的滤波结果如图6(b)、图6(c)和图6(d)所示。当r=3时，图6(b)中滤波结果包含有较多的噪声分量，当奇异值数量为r=2时和r=1时，信号中噪声分量得到90%的有效滤除。
综合上述仿真结果及分析表明，基于广义S变换和奇异值分解的心磁信号噪声自适应滤波消除设计方法是可行的。 

Claims (4)

1.心磁信号噪声自适应滤波消除设计方法，其特征在于，包括以下步骤： (1) 对实时采集到的心磁信号进行S变换，获得信号时频域内特征矩阵S变换时频矩阵； (2) 根据S变换时频矩阵时频分辨率的实时要求高低，在S变换中引入时频调节因子                                               
Figure 2013101425659100001DEST_PATH_IMAGE002
、
Figure 2013101425659100001DEST_PATH_IMAGE004
构建广义S变换，利用广义S变换时频调节因子调节信号时频特征的时频分辨率；(3) 采用奇异值分解方法分解广义S变换二维时频特征矩阵，得到左特征时域矩阵、右特征频域矩阵及特征值组成的对角阵，将单个奇异值对角阵分别与左右特征时域矩阵相乘得到单奇异值S变换时频矩阵，单奇异值S变换时频矩阵经反S变换即可得到单个主要奇异值所对应的信号时域成份；(4) 计算单个奇异值时域信号与噪声的相关系数，取相关系数r<0.5且单奇异值占总奇异值比例大于10%的奇异值所对应奇异矩阵区域为有效心磁信号奇异值分布区域，对有效心磁信号时频域进行反S变换重构信号从而实现心磁信号自适应时频滤波。
2.根据权利要求1所述的心磁信号噪声自适应滤波消除设计方法，其特征在于，所述步骤（1）中，所述心磁信号是指工作在无磁屏蔽室条件下利用高温射频超导量子干涉仪系统所测量到的心脏电心理过程在躯干上面产生的时变磁场强度信号；
所述步骤（1）中，所述S变换是指以高斯函数为窗函数的信号短时傅里叶变换，其表达式为
Figure 2013101425659100001DEST_PATH_IMAGE006
，式中，
Figure 2013101425659100001DEST_PATH_IMAGE008
为时间，
Figure 2013101425659100001DEST_PATH_IMAGE010
为频率，
Figure 2013101425659100001DEST_PATH_IMAGE012
为时移因子，
Figure 2013101425659100001DEST_PATH_IMAGE014
为虚数，
Figure 2013101425659100001DEST_PATH_IMAGE016
为表示实时连续心磁信号，
Figure 2013101425659100001DEST_PATH_IMAGE018
为高斯窗函数，其中
Figure 2013101425659100001DEST_PATH_IMAGE020
。
3.根据权利要求2所述的心磁信号噪声自适应滤波消除设计方法，其特征在于，所述步骤（2）中，所述广义S变换是指在S变换中引入时频调节因子、
Figure 902950DEST_PATH_IMAGE004
，其表达式为
Figure 2013101425659100001DEST_PATH_IMAGE024
；其中，
Figure 131675DEST_PATH_IMAGE002
为高斯窗幅度拉伸因子，为频率尺度拉伸因子，
Figure 244305DEST_PATH_IMAGE008
为时间，
Figure 681846DEST_PATH_IMAGE010
为频率，为时移因子，
Figure 646708DEST_PATH_IMAGE014
为虚数，为表示实时连续心磁信号。
4.根据权利要求3所述的心磁信号噪声自适应滤波消除设计方法，其特征在于，所述步骤（3）中，所述奇异值分解方法是指一种非退化正交矩阵分解法，具体是指心磁信号经广义S变换后的二维时频特征矩阵正交分解方法。
CN201310142565.9A 2013-04-23 2013-04-23 心磁信号噪声自适应滤波消除设计方法 Active CN103190898B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201310142565.9A CN103190898B (zh) 2013-04-23 2013-04-23 心磁信号噪声自适应滤波消除设计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201310142565.9A CN103190898B (zh) 2013-04-23 2013-04-23 心磁信号噪声自适应滤波消除设计方法

Publications (2)

Publication Number Publication Date
CN103190898A true CN103190898A (zh) 2013-07-10
CN103190898B CN103190898B (zh) 2014-09-10

Family

ID=48713848

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201310142565.9A Active CN103190898B (zh) 2013-04-23 2013-04-23 心磁信号噪声自适应滤波消除设计方法

Country Status (1)

Country Link
CN (1) CN103190898B (zh)

Cited By (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105527618A (zh) * 2016-02-26 2016-04-27 中国矿业大学(北京) 一种探地雷达地埋目标有效信号增强方法
CN108151873A (zh) * 2017-12-26 2018-06-12 广东石油化工学院 一种分离离心泵振动信号和循环水换热器振动信号的方法
CN109541455A (zh) * 2018-12-03 2019-03-29 国网江苏省电力有限公司南京供电分公司 一种基于s变换时频谱svd降噪的oltc冲击特性提取方法
CN110018429A (zh) * 2019-03-29 2019-07-16 中国科学院电子学研究所 一种消除磁探测平台振动干扰磁场的方法和系统
CN110646203A (zh) * 2019-08-23 2020-01-03 中国地质大学（武汉） 基于奇异值分解和自编码器的轴承故障特征提取方法
CN111012317A (zh) * 2020-01-18 2020-04-17 四川知周光声医疗科技有限公司 一种光声乳腺的图像重建方法及系统
CN111345806A (zh) * 2018-12-21 2020-06-30 四川锦江电子科技有限公司 一种心动周期检测方法及装置
CN112405072A (zh) * 2020-11-11 2021-02-26 上海交通大学 一种机床切削颤振的在线监测方法及装置
CN113017640A (zh) * 2021-02-23 2021-06-25 安徽省立医院(中国科学技术大学附属第一医院) 一种心磁信号背景噪声s变换域去除方法及系统
CN113180680A (zh) * 2021-05-17 2021-07-30 复旦大学 一种改进的基于奇异谱分析的心电信号降噪方法
CN113317793A (zh) * 2021-06-11 2021-08-31 宁波大学 心磁高频信号分析方法、存储介质及电子设备
CN113765502A (zh) * 2021-08-19 2021-12-07 宁波力斗智能技术有限公司 一种基于s域紧致奇异值分解的pd源滤波方法
CN115018018A (zh) * 2022-08-05 2022-09-06 北京航空航天大学杭州创新研究院 一种抑制背景噪声的双重空间滤波方法
CN116562202A (zh) * 2023-07-11 2023-08-08 广汽埃安新能源汽车股份有限公司 一种滤波组件分析方法及装置
CN117100276A (zh) * 2023-10-23 2023-11-24 山东大学齐鲁医院 心功能检测系统、计算机存储介质及终端

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101393248A (zh) * 2008-06-18 2009-03-25 昆明理工大学 一种基于s变换的输电线路故障行波波头精确定位方法
CN102685053A (zh) * 2012-05-15 2012-09-19 北京航空航天大学 一种基于广义s变换的通信信号调制识别方法
US20130079622A1 (en) * 2011-01-31 2013-03-28 Chenyu Wu Denoise MCG Measurements

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101393248A (zh) * 2008-06-18 2009-03-25 昆明理工大学 一种基于s变换的输电线路故障行波波头精确定位方法
US20130079622A1 (en) * 2011-01-31 2013-03-28 Chenyu Wu Denoise MCG Measurements
CN102685053A (zh) * 2012-05-15 2012-09-19 北京航空航天大学 一种基于广义s变换的通信信号调制识别方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
赵莉等: "互补型自适应滤波器在心磁信号处理中的应用", 《物理学报》 *
赵莉等: "小波变换在心磁信号处理中的应用", 《物理学报》 *

Cited By (23)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105527618A (zh) * 2016-02-26 2016-04-27 中国矿业大学(北京) 一种探地雷达地埋目标有效信号增强方法
CN108151873A (zh) * 2017-12-26 2018-06-12 广东石油化工学院 一种分离离心泵振动信号和循环水换热器振动信号的方法
CN109541455A (zh) * 2018-12-03 2019-03-29 国网江苏省电力有限公司南京供电分公司 一种基于s变换时频谱svd降噪的oltc冲击特性提取方法
CN111345806B (zh) * 2018-12-21 2023-01-31 四川锦江电子医疗器械科技股份有限公司 一种心动周期检测方法及装置
CN111345806A (zh) * 2018-12-21 2020-06-30 四川锦江电子科技有限公司 一种心动周期检测方法及装置
CN110018429A (zh) * 2019-03-29 2019-07-16 中国科学院电子学研究所 一种消除磁探测平台振动干扰磁场的方法和系统
CN110018429B (zh) * 2019-03-29 2021-01-15 中国科学院电子学研究所 一种消除磁探测平台振动干扰磁场的方法和系统
CN110646203A (zh) * 2019-08-23 2020-01-03 中国地质大学（武汉） 基于奇异值分解和自编码器的轴承故障特征提取方法
CN111012317B (zh) * 2020-01-18 2022-10-25 中川新迈科技有限公司 一种光声乳腺的图像重建方法及系统
CN111012317A (zh) * 2020-01-18 2020-04-17 四川知周光声医疗科技有限公司 一种光声乳腺的图像重建方法及系统
CN112405072A (zh) * 2020-11-11 2021-02-26 上海交通大学 一种机床切削颤振的在线监测方法及装置
CN113017640A (zh) * 2021-02-23 2021-06-25 安徽省立医院(中国科学技术大学附属第一医院) 一种心磁信号背景噪声s变换域去除方法及系统
CN113017640B (zh) * 2021-02-23 2024-06-04 安徽省立医院(中国科学技术大学附属第一医院) 一种心磁信号背景噪声s变换域去除方法及系统
CN113180680A (zh) * 2021-05-17 2021-07-30 复旦大学 一种改进的基于奇异谱分析的心电信号降噪方法
CN113180680B (zh) * 2021-05-17 2022-06-21 复旦大学 一种改进的基于奇异谱分析的心电信号降噪方法
CN113317793A (zh) * 2021-06-11 2021-08-31 宁波大学 心磁高频信号分析方法、存储介质及电子设备
CN113765502A (zh) * 2021-08-19 2021-12-07 宁波力斗智能技术有限公司 一种基于s域紧致奇异值分解的pd源滤波方法
CN113765502B (zh) * 2021-08-19 2023-10-20 宁波力斗智能技术有限公司 一种基于s域紧致奇异值分解的pd源滤波方法
CN115018018A (zh) * 2022-08-05 2022-09-06 北京航空航天大学杭州创新研究院 一种抑制背景噪声的双重空间滤波方法
CN116562202A (zh) * 2023-07-11 2023-08-08 广汽埃安新能源汽车股份有限公司 一种滤波组件分析方法及装置
CN116562202B (zh) * 2023-07-11 2023-09-08 广汽埃安新能源汽车股份有限公司 一种滤波组件分析方法及装置
CN117100276A (zh) * 2023-10-23 2023-11-24 山东大学齐鲁医院 心功能检测系统、计算机存储介质及终端
CN117100276B (zh) * 2023-10-23 2024-01-12 山东大学齐鲁医院 心功能检测系统、计算机存储介质及终端

Also Published As

Publication number Publication date
CN103190898B (zh) 2014-09-10

Similar Documents

Publication Publication Date Title
CN103190898B (zh) 心磁信号噪声自适应滤波消除设计方法
CN102697493B (zh) 一种快速的脑电信号中眼电伪迹自动识别和去除的方法
CN107192553B (zh) 基于盲源分离的齿轮箱复合故障诊断方法
CN102697495B (zh) 基于总体平均经验模式分解的二代小波肌电信号消噪方法
CN107144879B (zh) 一种基于自适应滤波与小波变换结合的地震波降噪方法
WO2006034024A2 (en) Method for adaptive complex wavelet based filtering of eeg signals
CN106771905B (zh) 一种适用于高频电流局部放电检测的脉冲提取方法
CN105974468B (zh) 一种能够同时进行五维地震数据重建和噪声压制的方法
CN104688220A (zh) 一种去除脑电信号中眼电伪迹的方法
CN109589114A (zh) 基于ceemd和区间阈值的肌电消噪方法
CN103961092A (zh) 基于自适应阈值处理的脑电信号去噪方法
CN105286860A (zh) 一种基于双树复小波能量差的运动想象脑电信号识别方法
CN114469124B (zh) 一种运动过程中异常心电信号的识别方法
CN104485113A (zh) 一种多故障源声发射信号分离方法
Figueiredo et al. Wavelet decomposition and singular spectrum analysis for electrical signal denoising
CN110161491A (zh) 一种针对微弱生命体的测距和呼吸频率估计方法
CN113640891A (zh) 一种基于奇异谱分析的瞬变电磁探测数据噪声滤除方法
Wang et al. Research on denoising algorithm for ECG signals
CN109567792A (zh) 一种单通道腹部记录胎儿心电图提取方法
Cao et al. A method for extracting weak impact signal in NPP based on adaptive Morlet wavelet transform and kurtosis
Ma et al. A multichannel nonlinear adaptive noise canceller based on generalized FLANN for fetal ECG extraction
Bandi et al. Single-channel fetal ECG extraction method based on extended Kalman filtering with singular value decomposition algorithm
CN105640505A (zh) 一种基于ar模型谱估计的脉搏信号随机噪声去噪方法
CN104778342B (zh) 一种基于小波奇异熵的心音特征提取方法
CN104935292A (zh) 一种基于源数估计的表面肌电信号自适应滤波方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant
TR01 Transfer of patent right

Effective date of registration: 20210924

Address after: Haiqi street, Ningbo City, Zhejiang Province

Patentee after: Ningbo Lidou Intelligent Technology Co.,Ltd.

Address before: 410082 School of electrical and information engineering, Hunan University, Changsha City, Hunan Province

Patentee before: He Yigang

TR01 Transfer of patent right