一种基于ar模型谱估计的脉搏信号随机噪声去噪方法
【技术领域】
[0001]本发明涉及一种基于AR模型谱估计的脉搏信号随机噪声去噪方法，属于信息科学 与医学融合技术领域。
【背景技术】
[0002] 人体内部各个生理系统之间是相互耦合的。反映人身体健康状态相对最重要、最 全面的是心脏血液循环系统，因此通过采集脉搏波进而分析心脏循环系统功能，能较全面 地反映人体的健康情况。但从人体采集到的脉搏波信号，由于脉搏波信号信噪比比较低，给 后续参数的准确测量带来了困难，所以对于噪声干扰的去除是非常重要而必须的。
[0003] 常用的脉搏信号去噪方法主要包括仅在时域内进行简单降噪处理的方法或传统 的建立在傅立叶变换基础上的去噪方法。仅在时域内进行简单降噪处理，即采用低通、带 通、高通等常规的滤波方法对脉搏波信号进行滤波去噪，然而常用的滤波器截止频率固定， 在噪声频率超过其截止频率时，无法消除噪声，截止频率过高会滤除有用信息。当噪声与信 号的频谱相近或重叠时，传统的单纯时域滤波或频域滤波往往无法达到很好的效果。传统 的建立在傅立叶变换基础上的去噪方法，在提高信噪比和提高空间分辨率两项指标上存在 矛盾。低通滤波能通过平滑抑制噪声，但同时会使信号的边沿变模糊。高通滤波可以使信号 边沿更加的陡峭，但背景噪声也同时被加强了。小波变换作为一种新型的时域分析法，由于 其具有良好的时频局部性，并且有快速算法(Mallat算法)加以实现，因而在去噪领域受到 了越来越多的这类方法处理比较简单，运算量较小，基于小波变换的信号去噪方法，一般有 阈值法、空域相关滤波法和模极大值法。但采用阈值法对脉搏波信号进行去噪，当信号具有 突变的不连续点或信噪比较低时，阈值法去噪会出现伪吉布斯现象，即在不连续点附近的 信号会在原有信号电平上上下跳变。空域相关滤波它能够保留信号突变点位置携带的重要 信息，从而对噪声有很好的滤波效果，但是传统的空域相关滤波算法很容易引入随机噪声。 基于模极大值的小波去噪算法，小波分解尺度的选择非常重要，尺度过小，小波系数受噪声 影响非常大，产生许多伪极值点，尺度过大又会使信号丢失某些重要的局部奇异性，因此需 要选择合适的尺度。
【发明内容】
[0004] 本发明提供了一种基于AR模型谱估计的脉搏信号随机噪声去噪方法，通过估计脉 搏信号功率信噪比，进而建立随机噪声信号模型，提出了利用AR模型谱估计，在频域内去除 随机噪声的方法，用于解决有用信号和噪声信号在频域宽度内重合时，用传统时域滤波和 频域滤波方法无法较好分离有用信号和噪声信号的问题;也用于解决基于小波变换的去噪 方法会产生伪吉布斯现象或引入随机噪声的问题;且本发明中，确定AR参数模型的阶数为 未补零前的信号长度，参数的选择简单且固定，但去噪效果依然很好，解决了参数选择困难 及因参数选择不当而影响到去噪效果的问题。
[0005] 本发明的技术方案是:一种基于AR模型谱估计的脉搏信号随机噪声去噪方法，所 述基于AR模型谱估计的脉搏信号随机噪声去噪方法的具体步骤如下：
[0006] Stepl、将采集到的脉搏信号分为长度相同的两段来计算信号的功率信噪比，通过 计算得到的信号功率信噪比估计噪声方差，进而建立随机噪声信号的模型，设计随机噪声 信号；
[0007] Step2、分别给两段信号补零使两段信号的大小均为最接近原信号长度的2的整数 次方，对补零后的混合信号进行傅里叶变换保留其相位谱；
[0008] Step3、对补零后的两段混合信号建立AR模型，确定AR模型及模型阶数，根据模型 求出模型的参数，代入功率谱密度估计公式估计出补零后的混合信号的功率谱;其中模型 阶数选择未补零前的每段信号的长度；
[0009] Step4、对所设计的随机噪声信号补零，使噪声信号的长度与补零后的混合信号的 长度相同，用与步骤Step3相同的方法估计出补零后随机噪声信号的功率谱，模型阶数也与 步骤Step3中相同；
[0010] Step5、用补零后的混合信号的功率谱减去补零后随机噪声信号的功率谱，得到有 效信号的功率谱，结合去噪前混合信号的相位谱，通过变换得到时域脉搏有效信号。
[0011] 所述步骤Stepl的具体步骤如下：
[0012] Stepl · 1、将脉搏信号x(m)(m=l，2,3，···，M)分为长度相同的两段信号xi(i)和X2
，使得X
且记/ = Au
[0013] Stepl .2、采集到的脉搏信号X(m) (m=l，2,3,'"，M)是由脉搏有效信号S(HI)和噪 声信号11(111)组成的，即叉(111)=8(111)+11(111)，贝1^1(；〇 = 81(；〇+111(；〇，叉2(：〇 = 82(：〇+112(：〇，取叉1 (i)、x2(i)的傅里叶变换，得到XKk)和x2(k)(k=0，l，…，I)，然后分别取它们幅值的平方， 并除以信号的长度I，作为对X i ( i )和X 2 ( i )真实功率谱的估计，即
1则信号X i ( i )、X 2 ( i )的功率分别为
则每段混合信号的平均功率为€ +A2)， 通过公式⑷二4认#)_<(/、')+^)，.\':(々)），估计出信号幻〇)、幻（〇的互功率谱， 因为X1(i)和X2(i)的互功率谱就是有效信号的功率谱，则脉搏有效信号的平均功率
KA·))，因为X1(i)、X2(i)这两段信号分别的总功率谱 就是有效信号与随机噪声信号功率谱之和，故可得噪声功率为G ，可求出信号的 功率信噪比为S/V/i二;
[0014] 其中，由于X1(i)、X2(i)之间有效信号相关性强、有效信号与随机噪声信号不相关、 随机噪声信号之间也不相关，故每段信号的自相关函数Λ (「）（i = l，2)为脉搏有效信号的 自相关函数rs(W与噪声自相关函数Γη(τ)之和，而两段信号的互相关函数(r)等于脉搏 有效信号的自相关函数rs(τ)，即/乂⑴=rs⑴+ r" (τ)，rv.2⑴=rs (r)，/ = 1，2，因此，X1 (i)、 x2(i)这两段信号分别的总功率谱Px(k)就是有效信号Ps(k)与噪声功率谱P n(k)之和，而它 们之间的互功率谱(幻就是有效信号的功率谱Ps(k)，即Px(k)=P s(k)+Pn(k)， Pvi(A)=P(A)-
[0015] 3七6?1.3、在时域中，脉搏信号11(1)32(1)的功率为
故每段信号的平均功率夕
[0016] Stepl .4、选择均值为0、方差为〇2、长度为I的服从高斯分布的伪白噪声序列作为 随机噪声信号的模型；由于白噪声的功率是用其方差来定义的，即Gu = 〇2，Gu表示白噪声的 功率，故采集到的脉搏信号中所含的随机噪声的功率为Gs_n = Gu = σ2，因此，SNR = Gs_s/Gs_n = (Gs_x-Gs_n)/Gs_ n=(Gs_x-Gu)/Gui% = G :,/d+l>，这样就可以求得白噪声模型的方差 σ2,进而建立出均值为0、方差为σ2、长度为I的伪高斯白噪声模型n(i) 了。
[0017] 所述步骤Step2的具体步骤如下：
[0018] 分别对脉搏信号"!^(!!!二^一^-^丨分成的两段长度相同的信号^^八^^)# 零得到新的信号幻"）、^")进行傅里叶变换，得到XKk)和X 2(k)(k = 0，l，…，N-1)，保留其 相位谱队(幻，其相位谱为傅里叶变换结果的虚部与实部之商，即为
[0019] (i) = imcig(X\{k))IrealiXi(k)), φΧι (k) = imag{X2(ky}/rml(X2(k))；：
[0020] 其中，n = l，2,3,…，N，N为最接近原信号x(m)长度M的2的整数次方，且N满足N 2 M。
[0021 ] 所述步骤Step3的具体步骤如下：
[0022] Step3.1、对脉搏信号X(m)(m=l，2,3，'"，M)分成的两段长度相同信号的其中一段 信号X1(i)补零后得到的信号X1(n)(n = l，2,3,…，N)，建立AR模型，N为最接近原信号x(m)长 度M的2的整数次方，且N满足N 2 M，建立的AR模型可由公另
表达，确定AR参数模型的阶数p为未补零前信号^(1)的长度，即p = I;
[0023] Step3.2、估计出信号χι(η)的前p + 1个自相关函数 /'(〇)，·(\(!)，弋(2)，…，
[0024] Step3.3、建立AR模型的参数ak和χΚη)的前p+1个自相关函数的关系，即建立AR模 型的正则方程，又称为Yu I e-Wa I ker方程，建立的方程如下：
[0026] 用心(0)，M1)，；V2)，___，^ (夕)代替方程中的《(〇)，~⑴，·_'，iV(A)，用 Levinson-Durbin递推算法求解方程，求出真实参数的估计值，即士，?2,…，<，及'：
[0027] Step3.4、将估计值\<，.沒「代入公式
中， 得到Χ1 (η)的功率谱尽(，）的估计，即
，对ω在单位圆上 均勾抽样，取分点个数为Xi(η)长度的4倍，记分点个数为nfft, 即nfft = 4*Ν，则得到离散谱
，式中：
=0,1，…，nfft_l;
[0028] Step3.5、对X2(n)进行步骤Step3.1至步骤Step3.4的同样处理，求出脉搏信号X (111)(111=1，2,3，一，)分成的两段长度相同信号的其中一段信号12(1)补零后得到的信号1 2 M- M (η )的功率谱估计☆ gV、，其中，xi(i)和 X2(i) (i = 1，2,3,'··,[τ])，. XI ( i ) = X ( m )
[0029] 所述步骤Step3.2的具体步骤如下：
[0030] Step3.2.1 JixKn)补零使其长度为最接近2N-1的2的整数次方，记补零后的数据 长度为N1，且N1应满足N 1 2 2N-1，补零后的信号记为λ·ινι(?);
[0031] Step3 · 2 · 2、对.?,(力)作傅里叶变换得尤以|(幻（^〇,1，...，1-1):，求11^ 1(幻的幅平 方，然后除以N1，得
[0032] Step3.2.3、对尤叫(々）作傅里叶逆变换，得xLVl(?)自相关函数的估计 而，⑴并不简单地等于自相关函数的估计弋(0，此时;\的为(Γ)的前N项与后N-I项的 组合，且后N-I项在前，前N项在后，即弋(r) = [& (rXM - iV + 2 : JV1);
[0033] ~ (r)(l: #)] = [CW1 - 況 + 2 : TV1):C(1: JV)]，τ = -(N-I)，-(N-2)…，-1，0，1，…，N-2， N-I;由此，得到& (〇: /?)=弋(〇 : /3) = c(l: P + 1)，因而估计出了信号Xl(n)的