技术领域
[0001] 本
发明涉及一种基于AR模型谱估计的脉搏信号随机噪声去噪方法，属于信息科学与医学融合技术领域。
背景技术
[0002] 人体内部各个生理系统之间是相互耦合的。反映人身体
健康状态相对最重要、最全面的是心脏血液循环系统，因此通过采集脉搏波进而分析心脏循环系统功能，能较全面地反映人体的健康情况。但从人体采集到的脉搏波信号，由于脉搏波信号
信噪比比较低，给后续参数的准确测量带来了困难，所以对于噪声干扰的去除是非常重要而必须的。
[0003] 常用的脉搏信号去噪方法主要包括仅在时域内进行简单降噪处理的方法或传统的建立在
傅立叶变换基础上的去噪方法。仅在时域内进行简单降噪处理，即采用低通、带通、高通等常规的滤波方法对脉搏波信号进行滤波去噪，然而常用的
滤波器截止
频率固定，在噪声频率超过其截止频率时，无法消除噪声，截止频率过高会滤除有用信息。当噪声与信号的
频谱相近或重叠时，传统的单纯时域滤波或频域滤波往往无法达到很好的效果。传统的建立在傅立叶变换基础上的去噪方法，在提高信噪比和提高空间
分辨率两项指标上存在矛盾。低通滤波能通过平滑抑制噪声，但同时会使信号的边沿变模糊。高通滤波可以使信号边沿更加的陡峭，但背景噪声也同时被加强了。
小波变换作为一种新型的时域分析法，由于其具有良好的时频局部性，并且有快速
算法(Mallat算法)加以实现，因而在去噪领域受到了越来越多的这类方法处理比较简单，运算量较小，基于小波变换的信号去噪方法，一般有
阈值法、
空域相关滤波法和模极大值法。但采用阈值法对脉搏波信号进行去噪，当信号具有突变的不连续点或信噪比较低时，阈值法去噪会出现伪吉布斯现象，即在不连续点附近的信号会在原有信号电平上上下跳变。空域相关滤波它能够保留信号突变点
位置携带的重要信息，从而对噪声有很好的滤波效果，但是传统的空域相关滤波算法很容易引入随机噪声。基于模极大值的小波去噪算法，小波分解尺度的选择非常重要，尺度过小，小波系数受噪声影响非常大，产生许多伪极值点，尺度过大又会使信号丢失某些重要的局部奇异性，因此需要选择合适的尺度。
发明内容
[0004] 本发明提供了一种基于AR模型谱估计的脉搏信号随机噪声去噪方法，通过估计脉搏信号功率信噪比，进而建立随机噪声信号模型，提出了利用AR模型谱估计，在频域内去除随机噪声的方法，用于解决有用信号和噪声信号在频域宽度内重合时，用传统时域滤波和频域滤波方法无法较好分离有用信号和噪声信号的问题；也用于解决基于小波变换的去噪方法会产生伪吉布斯现象或引入随机噪声的问题；且本发明中，确定AR参数模型的阶数为未补零前的信号长度，参数的选择简单且固定，但去噪效果依然很好，解决了参数选择困难及因参数选择不当而影响到去噪效果的问题。
[0005] 本发明的技术方案是：一种基于AR模型谱估计的脉搏信号随机噪声去噪方法，所述基于AR模型谱估计的脉搏信号随机噪声去噪方法的具体步骤如下：
[0006] Step1、将采集到的脉搏信号分为长度相同的两段来计算信号的功率信噪比，通过计算得到的信号功率信噪比估计噪声方差，进而建立随机噪声信号的模型，设计随机噪声信号；
[0007] Step2、分别给两段信号补零使两段信号的大小均为最接近原信号长度的2的整数次方，对补零后的混合信号进行傅里叶变换保留其
相位谱；
[0008] Step3、对补零后的两段混合信号建立AR模型，确定AR模型及模型阶数，根据模型求出模型的参数，代入
功率谱密度估计公式估计出补零后的混合信号的功率谱；其中模型阶数选择未补零前的每段信号的长度；
[0009] Step4、对所设计的随机噪声信号补零，使噪声信号的长度与补零后的混合信号的长度相同，用与步骤Step3相同的方法估计出补零后随机噪声信号的功率谱，模型阶数也与步骤Step3中相同；
[0010] Step5、用补零后的混合信号的功率谱减去补零后随机噪声信号的功率谱，得到有效信号的功率谱，结合去噪前混合信号的相位谱，通过变换得到时域脉搏有效信号。
[0011] 所述步骤Step1的具体步骤如下：
[0012] Step1.1、将脉搏信号x(m)(m＝1,2,3,…,M)分为长度相同的两段信号x1(i)和x2(i)， 使得x1(i)＝x(m)， x2(i)＝x(m)，且记
[0013] Step1.2、采集到的脉搏信号x(m)(m＝1,2,3,…,M)是由脉搏有效信号s(m)和噪声信号n(m)组成的，即x(m)＝s(m)+n(m)，则x1(i)＝s1(i)+n1(i)，x2(i)＝s2(i)+n2(i)，取x1(i)、x2(i)的傅里叶变换，得到X1(k)和X2(k)(k＝0,1,…,I)，然后分别取它们幅值的平方，并除以 信号的 长度I ，作为对x 1( i )和x2( i )真实功率 谱的 估计 ，即则信号x1( i )、x2(i )的功率分别为
则每段混合信号的
平均功率为
通过公式 估计出信号x1(i)、x2(i)的互功率谱，
因为x1(i)和x2(i)的互功率谱就是有效信号的功率谱，则脉搏有效信号的平均功率因为x1(i)、x2(i)这两段信号分别的总功率谱
就是有效信号与随机噪声信号功率谱之和，故可得噪声功率为 可求出信号的
功率信噪比为
[0014] 其中，由于x1(i)、x2(i)之间有效信号相关性强、有效信号与随机噪声信号不相关、随机噪声信号之间也不相关，故每段信号的自相关函数 (i＝1,2)为脉搏有效信号的自相关函数rs(τ)与噪声自相关函数rn(τ)之和，而两段信号的互相关函数 等于脉搏有效信号的自相关函数rs(τ)，即 因此，x1(i)、x2(i)这两段信号分别的总功率谱Px(k)就是有效信号Ps(k)与噪声功率谱Pn(k)之和，而它们之间的互功率谱 就是有效信号的功率谱Ps(k)，即Px(k)＝Ps(k)+Pn(k)，
[0015] Step1.3、在时域中，脉搏信号x1(i)、x2(i)的功率为故每段信号的平均功率为
[0016] Step1.4、选择均值为0、方差为σ2、长度为I的服从高斯分布的伪白噪声序列作为随机噪声信号的模型；由于白噪声的功率是用其方差来定义的，即Gu＝σ2，Gu表示白噪声的2
功率，故采集到的脉搏信号中所含的随机噪声的功率为Gs_n＝Gu＝σ，因此，SNR＝Gs_s/Gs_n＝(Gs_x-Gs_n)/Gs_n＝(Gs_x-Gu)/Gu，故 这样就可以求得白噪声模型的方差
σ2，进而建立出均值为0、方差为σ2、长度为I的伪高斯白噪声模型n(i)了。
[0017] 所述步骤Step2的具体步骤如下：
[0018] 分别对脉搏信号x(m)(m＝1,2,3,…,M)分成的两段长度相同的信号x1(i)、x2(i)补零得到新的信号x1(n)、x2(n)进行傅里叶变换，得到X1(k)和X2(k)(k＝0,1,…,N-1)，保留其相位谱 其相位谱为傅里叶变换结果的
虚部与
实部之商，即为
[0019]
[0020] 其中，n＝1,2,3,…，N，N为最接近原信号x(m)长度M的2的整数次方，且N满足N≥M。
[0021] 所述步骤Step3的具体步骤如下：
[0022] Step3.1、对脉搏信号x(m)(m＝1,2,3,…,M)分成的两段长度相同信号的其中一段信号x1(i)补零后得到的信号x1(n)(n＝1,2,3,…,N)，建立AR模型，N为最接近原信号x(m)长度M的2的整数次方，且N满足N≥M，建立的AR模型可由公式表达，确定AR参数模型的阶数p为未补零前信号x1(i)的长度，即p＝I；
[0023] S t e p 3 . 2 、估 计 出 信 号 x 1 ( n )的 前 p + 1 个 自 相 关 函 数[0024] Step3.3、建立AR模型的参数ak和x1(n)的前p+1个自相关函数的关系，即建立AR模型的正则方程，又称为Yule-Walker方程，建立的方程如下：
[0025]
[0026] 用 代替方程中的 用Levinson-Durbin递推算法求解方程，求出真实参数的估计值，即
[0027] Step3.4、将估计值 代入公式 中，得到x1(n)的功率谱 的估计，即 对ω在单位圆上
均匀抽样，取分点个数为x1(n)长度的4倍，记分点个数为nfft,即nfft＝4*N，则得到离散谱式中， 而 l
＝0,1,…,nfft-1；
[0028] Step3.5、对x2(n)进行步骤Step3.1至步骤Step3.4的同样处理，求出脉搏信号x(m)(m＝1,2,3,…,M)分成的两段长度相同信号的其中一段信号x2(i)补零后得到的信号x2(n)的功率谱估计 其中，x1(i)和x2(i) x1(i)＝x(m)x2(i)＝x(m) 且记
[0029] 所述步骤Step3.2的具体步骤如下：
[0030] Step3.2.1、对x1(n)补零使其长度为最接近2N-1的2的整数次方，记补零后的数据长度为N1，且N1应满足N1≥2N-1，补零后的信号记为
[0031] Step3.2.2、对 作傅里叶变换得 求 的幅平方，然后除以N1，得
[0032] Step3.2.3、对 作傅里叶逆变换，得 自相关函数的估计而 并不简单地等于自相关函数的估计 此时 为 的前N项与后N-1项的
组合，且后N-1项在前，前N项在后，即
[0033] τ＝-(N-1),-(N-2)…,-1,0,1,…,N-2,N-1；由此，得到 因而估计出了信号x1(n)的前p+1个自相关
函数
[0034] 所述步骤Step4的具体步骤如下：
[0035] 对所建立的服从高斯分布的伪白噪声序列的随机噪声信号的模型，即伪高斯白噪声信号n(i)补零，得到信号n(n)，使噪声信号的长度与混合信号x1(n)、x2(n)的长度相同，即长度为N，n(n)，n＝1,2,3,…,N，用与步骤Step3中相同的方法，估计出噪声信号的功率谱模型阶数和分点个数也与步骤Step3中的选择相同；
[0036] 其中，混合信号x1(n)、x2(n)为脉搏信号x(m)(m＝1,2,3,…,M)分为长度相同的两段信号x1(i)和x2(i)， 补零得到的新的信号，n＝1,2,3,…,N，N为最接近原信号x(m)长度M的2的整数次方，且N满足N≥M。
[0037] 所述步骤Step5的具体步骤如下：
[0038] Step5.1、信号x1(n)的功率谱估计 减去信号n(n)的功率谱估计得到信号s1(n)的功率谱估计 即 将 记为
[0039] Step5.2、对 作逆变换，得 将 的前 项与后 项调换顺序，使后 项在前，前 项在后，得到新的序列 即
和真实自相关函数 有如下关系
[0040]
[0041] 求得 为 中间的2p+1项，即记 为
[0042] Step5.3、对 补零得到 使 的长度与x1(n)相同，为N1，且使 的前p+1项等于的前p+1项，即 使 的后p+1项等于 的后p+1项，即使 的中间N1-2( p+1 )项的值均等于0 ，即
取 然后对 进行傅里叶变换，则变换结果的实部为有
2
效信号s1(n)的傅里叶变换幅值的平方|S1(k)，对变换结果的实部开平方，得到有效信号s1(n)傅里叶变换的幅值|S1(k)|；
[0043] Step5.4、结合信号x1(n)的相位谱 得到信号s1(n)的傅里叶变换结果S1(k)，即 对S1(k)取傅里叶逆变换，得到有效信号s1(n)，且取信号的前I个值，即可得到包含在混合信号x1(i)中的有效信号s1(i)；
[0044] Step5.5、用与步骤Step5.1至Step5.4相同的方法可以得到s2(i),将s2(i)与s1(i)首尾相接，即可得到去噪后的整段信号s(m)；
[0045] 其中，信号x1(n)、x2(n)为脉搏信号x(m)(m＝1,2,3,…,M)分为长度相同的两段信号x1(i)和x2(i)， 补零得到的新的信号，n＝1,2,…,N，N为最接近原信号x(m)长度M的2的整数次方，且N满足N≥M。
[0046] 本发明的有益效果是：通过估计脉搏信号功率信噪比，进而建立随机噪声信号模型，提出了利用AR模型谱估计，在频域内去除随机噪声的方法，根据不同的信号类型及数据长度，确定混合信号和随机噪声信号的AR模型及模型阶数，然后分别构建AR模型，求出模型参数，再根据模型参数求出混合信号和随机噪声信号的功率谱估计，进而在频率域去除随机干扰，最后变换到时域得到去除随机噪声后的有用信号，去噪后的信号更加平滑，在明显去除了随机噪声的基础上极大程度地保留了原信号的形状；再则，谱估计是信号内在本质的一种表现形式，利用AR模型谱估计方法去除随机噪声后的信号，既不会出现伪吉布斯现象又不会引入随机噪声,同时也没有降
低信号的分辨率及保真度；最后，本方法中，参数的选择主要涉及AR模型阶数的选择，可以确定AR模型的阶数为未补零前的信号长度，对于任意不同长度的待处理信号皆可做同样的选择，且不影响谱估计的效果，使得该方法简便易行。
附图说明
[0048] 图2为本发明中计算功率信噪比SNR的流程图；
[0049] 图3为本发明中对脉搏信号x(m)(m＝1,2,3,…,M)分成的两段长度相同信号的其中一段信号x1(i)利用AR模型进行去噪的流程图；对另一段信号x2(i)进行去噪的流程相同；
[0050] 图4为所采集的一段女性的真实脉搏信号x(m)(m＝1,2,3,…,M)；脉搏数据长度M＝20000；
[0051] 图5为对图4所示脉搏信号x(m)(m＝1,2,3,…,M)去噪后的信号s(m)(m＝1,2,3,…,M)；数据长度M＝20000；
[0052] 图6为将图4所示真实脉搏信号x(m)(m＝1,2,3,…,M)分成的两段长度相同信号中的前半段信号x1(i)(i＝1,2,3,…,I)；脉搏数据长度I＝10000；
[0053] 图7为对图6所示的脉搏信号x1(i)(i＝1,2,3,…,I)进行去噪后的信号s1(i)(i＝1,2,3,…,I)；脉搏数据长度I＝10000；
[0054] 图8为将图4所示真实脉搏信号x(m)(m＝1,2,3,…,M)分成的两段长度相同信号中的后半段信号x2(i)(i＝1,2,3,…,I)；脉搏数据长度I＝10000；
[0055] 图9为对图8所示的脉搏信号x2(i)(i＝1,2,3,…,I)进行去噪后的信号s2(i)(i＝1,2,3,…,I)；脉搏数据长度I＝10000；
[0056] 图10为所建立的伪高斯白噪声信号n(i)(i＝1,2,3,…,I)；信号长度为I＝10000；该伪高斯白噪声信号的均值为0、方差σ2为755.2006；
[0057] 图11为信号x1(n)(n＝1,2,3,…,N)的离散功率谱估计 其中，x1(n)为对图6所示的脉搏信号x1(i)(i＝1,2,3,…,I)补零后的信号，N＝32768；
[0058] 图12为信号x2(n)(n＝1,2,3,…,N)的离散功率谱估计 其中，x2(n)为对图8所示的脉搏信号x2(i)(i＝1,2,3,…,I)补零后的信号，N＝32768；
[0059] 图13为信号n(n)(n＝1,2,3,…,N)的功率谱估计 其中，n(n)为对图10所示的伪高斯白噪声信号n(i)(i＝1,2,3,…,I)补零后的信号，N＝32768；
[0060] 图14为信号s1(n)(n＝1,2,3,…,N)的功率谱估计 其中，s1(n)为信号x1(n)(n＝1,2,3,…,N)去噪后的信号，N＝32768；
[0061] 图15为信号s1(n)(n＝1,2,3,…,N)的自相关函数 (τ＝1,2,…,K)；其中，s1(n)为信号x1(n)(n＝1,2,3,…,N)去噪后的信号，N＝32768，K＝65536；
[0062] 图16为信号s2(n)(n＝1,2,3,…,N)的功率谱估计 其中，s2(n)为信号x2(n)(n＝1,2,3,…,N)去噪后的信号，N＝32768；
[0063] 图17为信号s2(n)(n＝1,2,3,…,N)的自相关函数 (τ＝1,2,…,K)；其中，s2(n)为对信号x2(n)(n＝1,2,3,…,N)去噪后的信号，N＝32768，K＝65536；
[0064] 图18为在图6所示信号x1(i)(i＝1,2,3,…,I)中截取的一段信号x1_1(t)(t＝1,2,3,…,T)；其中，数据长度T＝3000，该段信号截取的是x1(i)中的前三千个
采样点数据；
[0065] 图19为在图7所示信号s1(i)(i＝1,2,3,…,I)中截取的一段信号s1_1(t)(t＝1,2,3,…,T)；其中，数据长度T＝3000，该段信号截取的是s1(i)中的前三千个采样点数据；
[0066] 图20为在图8所示信号x2(i)(i＝1,2,3,…,I)中截取的一段信号x2_1(t)(t＝1,2,3,…,T)；其中，数据长度T＝3000，该段信号截取的是x2(i)中的前三千个采样点数据；
[0067] 图21为在图9所示信号s2(i)(i＝1,2,3,…,I)中截取的一段信号s2_1(t)(t＝1,2,3,…,T)；其中，数据长度T＝3000，该段信号截取的是s2(i)中的前三千个采样点数据。
具体实施方式
[0068]
实施例1：如图1-21所示，一种基于AR模型谱估计的脉搏信号随机噪声去噪方法，本方法的步骤如下：
[0069] Step1、将采集到的脉搏信号分为长度相同的两段来计算信号的功率信噪比，通过计算得到的信号功率信噪比估计噪声方差，进而建立随机噪声信号的模型，设计随机噪声信号；
[0070] Step2、分别给两段信号补零使两段信号的大小均为最接近原信号长度的2的整数次方，对补零后的混合信号进行傅里叶变换保留其相位谱；
[0071] Step3、对补零后的两段混合信号建立AR模型，确定AR模型及模型阶数，根据模型求出模型的参数，代入功率谱密度估计公式估计出补零后的混合信号的功率谱；其中模型阶数选择未补零前的每段信号的长度；
[0072] Step4、对所设计的随机噪声信号补零，使噪声信号的长度与补零后的混合信号的长度相同，用与步骤Step3相同的方法估计出补零后随机噪声信号的功率谱，模型阶数也与步骤Step3中相同；
[0073] Step5、用补零后的混合信号的功率谱减去补零后随机噪声信号的功率谱，得到有效信号的功率谱，结合去噪前混合信号的相位谱，通过变换得到时域脉搏有效信号。
[0074] 所述步骤Step1的具体步骤如下：
[0075] Step1.1、将脉搏信号x(m)(m＝1,2,3,…,M)分为长度相同的两段信号x1(i)和x2(i)， 使得x1(i)＝x(m)， x2(i)＝x(m)，且记
[0076] Step1.2、采集到的脉搏信号x(m)(m＝1,2,3,…,M)是由脉搏有效信号s(m)和噪声信号n(m)组成的，即x(m)＝s(m)+n(m)，则x1(i)＝s1(i)+n1(i)，x2(i)＝s2(i)+n2(i)，取x1(i)、x2(i)的傅里叶变换，得到X1(k)和X2(k)(k＝0,1,…,I)，然后分别取它们幅值的平方，并除以 信号的 长度I ，作为对x 1( i )和x2( i )真实功率 谱的 估计 ，即则 信号x1( i )、x2( i )的功 率分别 为
则 每 段 混 合 信 号 的 平 均 功 率 为
通过公式 估计出信号x1
(i)、x2(i)的互功率谱，因为x1(i)和x2(i)的互功率谱就是有效信号的功率谱，则脉搏有效信号的平均功率 因为x1(i)、x2(i)这两段信号
分别的总功率谱就是有效信号与随机噪声信号功率谱之和，故可得噪声功率为
可求出信号的功率信噪比为
[0077] 其中，由于x1(i)、x2(i)之间有效信号相关性强、有效信号与随机噪声信号不相关、随机噪声信号之间也不相关，故每段信号的自相关函数 (i＝1,2)为脉搏有效信号的自相关函数rs(τ)与噪声自相关函数rn(τ)之和，而两段信号的互相关函数 等于脉搏有效信号的自相关函数rs(τ)，即 i＝1,2，因此，x1(i)、x2(i)这两段信号分别的总功率谱Px(k)就是有效信号Ps(k)与噪声功率谱Pn(k)之和，而它们之间的互功率谱 就是有效信号的功率谱Ps(k)，即Px(k)＝Ps(k)+Pn(k)，
[0078] Step1.3、在时域中，脉搏信号x1(i)、x2(i)的功率为故每段信号的平均功率为
[0079] Step1.4、选择均值为0、方差为σ2、长度为I的服从高斯分布的伪白噪声序列作为随机噪声信号的模型；由于白噪声的功率是用其方差来定义的，即Gu＝σ2，Gu表示白噪声的功率，故采集到的脉搏信号中所含的随机噪声的功率为Gs_n＝Gu＝σ2，因此，
[0080] SNR＝Gs_s/Gs_n＝(Gs_x-Gs_n)/Gs_n＝(Gs_x-Gu)/Gu，故 这样就可以2 2
求得白噪声模型的方差σ，进而建立出均值为0、方差为σ、长度为I的伪高斯白噪声模型n(i)了。
[0081] 所述步骤Step2的具体步骤如下：
[0082] 分别对脉搏信号x(m)(m＝1,2,3,…,M)分成的两段长度相同的信号x1(i)、x2(i)补零得到新的信号x1(n)、x2(n)进行傅里叶变换，得到X1(k)和X2(k)(k＝0,1,…,N-1)，保留其相位谱 其相位谱为傅里叶变换结果的虚部与实部之商，即为
[0083]
[0084] 其中，n＝1,2,3,…,N，N为最接近原信号x(m)长度M的2的整数次方，且N满足N≥M。
[0085] 所述步骤Step3的具体步骤如下：
[0086] Step3.1、对脉搏信号x(m)(m＝1,2,3,…,M)分成的两段长度相同信号的其中一段信号x1(i)补零后得到的信号x1(n)(n＝1,2,3,…,N)，建立AR模型，N为最接近原信号x(m)长度M的2的整数次方，且N满足N≥M，建立的AR模型可由公式表达，确定AR参数模型的阶数p为未补零前信号x1(i)的长度，即p＝I；
[0087] Step3.2、估计出信号x1(n)的前p+1个自相关函数
[0088] Step3.3、建立AR模型的参数ak和x1(n)的前p+1个自相关函数的关系，即建立AR模型的正则方程，又称为Yule-Walker方程，建立的方程如下：
[0089]
[0090] 用 代替方程中的 用Levinson-Durbin递推算法求解方程，求出真实参数的估计值，即
[0091] Step3.4、将估计值 代入公式 中，得到x1(n)的功率谱 的估计，即 对ω在单位圆上
均匀抽样，取分点个数为x1(n)长度的4倍，记分点个数为nfft,即nfft＝4*N，则得到离散谱式中， 而 l
＝0,1,…,nfft-1；
[0092] Step3.5、对x2(n)进行步骤Step3.1至步骤Step3.4的同样处理，求出脉搏信号x(m)(m＝1,2,3,…,M)分成的两段长度相同信号的其中一段信号x2(i)补零后得到的信号x2(n)的功率谱估计 其中，x1(i)和x2(i) x1(i)＝x(m)x2(i)＝x(m) 且记
[0093] 所述步骤Step3.2的具体步骤如下：
[0094] Step3.2.1、对x1(n)补零使其长度为最接近2N-1的2的整数次方，记补零后的数据长度为N1，且N1应满足N1≥2N-1，补零后的信号记为
[0095] Step3.2.2、对 作傅里叶变换得 (k＝0,1,…,N1-1)，求 的幅平方，然后除以N1，得
[0096] Step3.2.3、对 作傅里叶逆变换，得 自相关函数的估计而 并不简单地等于自相关函数的估计 此时 为 的前N项与后N-1项的
组合，且后N-1项在前，前N项在后，即
[0097] τ＝-(N-1),-(N-2)…,-1,0,1,…,N-2,N-1；由此，得到 因而估计出了信号x1(n)的前p+1个自相关
函数
[0098] 所述步骤Step4的具体步骤如下：
[0099] 对所建立的服从高斯分布的伪白噪声序列的随机噪声信号的模型，即伪高斯白噪声信号n(i)补零，得到信号n(n)，使噪声信号的长度与混合信号x1(n)、x2(n)的长度相同，即长度为N，n(n)，n＝1,2,3,…,N，用与步骤Step3中相同的方法，估计出噪声信号的功率谱模型阶数和分点个数也与步骤Step3中的选择相同；
[0100] 其中，混合信号x1(n)、x2(n)为脉搏信号x(m)(m＝1,2,3,…,M)分为长度相同的两段信号x1(i)和x2(i)， 补零得到的新的信号，n＝1,2,3,…,N，N为最接近原信号x(m)长度M的2的整数次方，且N满足N≥M。
[0101] 所述步骤Step5的具体步骤如下：
[0102] Step5.1、信号x1(n)的功率谱估计 减去信号n(n)的功率谱估计得到信号s1(n)的功率谱估计 即 将 记为
[0103] Step5.2、对 作逆变换，得 将 的前 项与后 项调换顺序，使后 项在前，前 项在后，得到新的序列 即
和真实自相关函数 有如下关
系
[0104]
[0105] 求得 为 中间的2p+1项，即记 为
[0106] Step5.3、对 补零得到 使 的长度与x1(n)相同，为N1，且使 的前p+1项等于的前p+1项，即 使 的后p+1项等于 的后p+1项，即使 的中间N1-2( p+1 )项的值均等于0 ，即
取 然后对 进行傅里叶变换，则变换结果的实部为有
效信号s1(n)的傅里叶变换幅值的平方|S1(k)2，对变换结果的实部开平方，得到有效信号s1(n)傅里叶变换的幅值|S1(k)|；
[0107] Step5.4、结合信号x1(n)的相位谱 得到信号s1(n)的傅里叶变换结果S1(k)，即 对S1(k)取傅里叶逆变换，得到有效信号s1(n)，且取信号的前I个值，即可得到包含在混合信号x1(i)中的有效信号s1(i)；
[0108] Step5.5、用与步骤Step5.1至Step5.4相同的方法可以得到s2(i),将s2(i)与s1(i)首尾相接，即可得到去噪后的整段信号s(m)；
[0109] 其中，信号x1(n)、x2(n)为脉搏信号x(m)(m＝1,2,3,…,M)分为长度相同的两段信号x1(i)和x2(i)， 补零得到的新的信号，n＝1,2,…,N，N为最接近原信号x(m)长度M的2的整数次方，且N满足N≥M。
[0110] 实施例2：如图1-21所示，一种基于AR模型谱估计的脉搏信号随机噪声去噪方法，本方法的步骤如下：
[0111] Step1、将采集到的脉搏信号分为长度相同的两段来计算信号的功率信噪比，通过计算得到的信号功率信噪比估计噪声方差，进而建立随机噪声信号的模型；
[0112] 所述步骤Step1的具体步骤如下：
[0113] Step1.1、采集到的一位女性的真实脉搏信号数据x(m)(m＝1,2,3,…,M)，如图4，脉搏数据长度M＝20000，将x(m)分为长度相同的两段信号x1(i)和x2(i)(i＝1,2,3,…I)，I＝10000，x1(i)和x2(i)分别如图6、8所示；
[0114] Step1.2、采集到的脉搏信号x(m)(m＝1,2,3,…,M)是由脉搏有效信号s(m)和噪声信号n(m)组成的，即x(m)＝s(m)+n(m)，则x1(i)＝s1(i)+n1(i)，x2(i)＝s2(i)+n2(i)，取x1(i)、x2(i)的傅里叶变换，得到X1(k)和X2(k)(k＝0,1,…,I)，然后分别取它们幅值的平方，并除以 信号的 长度I ，作为对x 1( i )和x2( i )真实功率 谱的 估计 ，即则信号x1( i )、x2( i )的功率分 别为
则 每 段 混 合 信 号 的 平 均 功 率 为
通过公式 估计出信号x1
(i)、x2(i)的互功率谱，因为x1(i)和x2(i)的互功率谱就是有效信号的功率谱，则脉搏有效信号的平均功率 因为x1(i)、x2(i)这两段信号
分别的总功率谱就是有效信号与随机噪声信号功率谱之和，故可得噪声功率为
可求出信号的功率信噪比为
[0115] 这里求得 SNR＝8.5035；
[0116] 其中，由于x1(i)、x2(i)之间有效信号相关性强、有效信号与随机噪声信号不相关、随机噪声信号之间也不相关，故每段信号的自相关函数 (i＝1,2)为脉搏有效信号的自相关函数rs(τ)与噪声自相关函数rn(τ)之和，而两段信号的互相关函数 等于脉搏有效信号的自相关函数rs(τ)，即 i＝1,2，因此，x1(i)、x2(i)这两段信号分别的总功率谱Px(k)就是有效信号Ps(k)与噪声功率谱Pn(k)之和，而它们之间的互功率谱 就是有效信号的功率谱Ps(k)，即Px(k)＝Ps(k)+Pn(k)，
[0117] Step1.3、在时域中，脉搏信号x1(i)、x2(i)的功率为故每段信号的平均功率为
[0118] Step1.4、选择均值为0、方差为σ2、长度为I的服从高斯分布的伪白噪声序列作为随机噪声信号的模型；则 求出σ2，这里求得σ2＝755.2006；由于白噪声2
的功率是用其方差来定义的，即Gu＝σ，Gu表示白噪声的功率，故采集到的脉搏信号中所含的随机噪声的功率为Gs_n＝Gu＝σ2，因此，SNR＝Gs_s/Gs_n＝(Gs_x-Gs_n)/Gs_n＝(Gs_x-Gu)/Gu，故这样就可以建立出均值为0、方差为σ2、长度为I的伪高斯白噪声模型n
(i)(i＝1,2,…,10000)了。建立的伪高斯白噪声信号n(i)，如图10所示；
[0119] Step2、分别给两段信号补零使两段信号的大小均为最接近原信号长度的2的整数次方，对补零后的混合信号进行傅里叶变换保留其相位谱；
[0120] 所述步骤Step2的具体步骤如下：
[0121] 分别对脉搏信号x(m)(m＝1,2,3,…,M)分成的两段长度相同的信号x1(i)、x2(i)补零得到新的信号x1(n)、x2(n)进行傅里叶变换，得到X1(k)和X2(k)(k＝0,1,…,N-1)，保留其相位谱 其相位谱为傅里叶变换结果的虚部与实部之商，即为
[0122]
[0123] 其中，n＝1,2,3，…,N，N为最接近原信号x(m)长度M的2的整数次方，且N满足N≥M，计算得N＝32768。
[0124] Step3、对补零后的两段混合信号建立AR模型，确定AR模型及模型阶数，根据模型求出模型的参数，代入功率谱密度估计公式估计出补零后的混合信号的功率谱；其中模型阶数选择未补零前的每段信号的长度；
[0125] 所述步骤Step3的具体步骤如下：
[0126] Step3.1、对脉搏信号x(m)(m＝1,2,3,…,M)分成的两段长度相同信号的其中一段信号x1(i)补零后得到的信号x1(n)(n＝1,2,3,…,N)，这里，N＝32768，建立AR模型，立的AR模型可由公式 表达，确定AR参数模型的阶数p为未补零前的信号的x1(i)长度，即p＝I＝10000；
[0127] Step3.2、估计出信号x1(n)的前p+1个自相关函数
[0128] Step3.3、建立AR模型的参数ak和x1(n)的前p+1个自相关函数的关系，即建立AR模型的正则方程，又称为Yule-Walker方程，建立的方程如下：
[0129]
[0130] 用 代替方程中的 用Levinson-Durbin递推算法求解方程，求出真实参数的估计值，即
[0131] Step3.4、将估计值 代入公式 中，得到x1(n)的功率谱 的估计，即 对ω在单位圆上均匀抽
样，取分点个数为长度的4倍，记分点个数为nfft,即nfft＝4*N＝131872，则得到离散谱[0132] 式中， 而
l＝0,1,…,nfft-1；x1(n)的功率谱估计 如图11所示；
[0133] Step3.5、对x2(n)进行步骤Step3.1至步骤Step3.4的同样处理，求出脉搏信号x(m)(m＝1,2,3,…,M)分成的两段长度相同信号的其中一段信号x2(i)补零后得到的信号x2(n)的功率谱估计 如图12所示；
[0134] Step4、对所建立的随机噪声信号的模型补零，使噪声信号的长度与补零后的混合信号的长度相同，用与步骤Step3相同的方法估计出噪声信号模型的功率谱，模型阶数也与步骤Step3中相同；
[0135] 所述步骤Step4的具体步骤如下：
[0136] 对所建立的服从高斯分布的伪白噪声序列的随机噪声信号即伪高斯白噪声信号n(i)补零，得到信号n(n)，使噪声信号的长度与混合信号x1(n)、x2(n)的长度相同，即长度为N，n(n)(n＝1,2,3,…,N)，用与步骤Step3中相同的方法，估计出噪声信号的功率谱模型阶数和分点个数也与步骤Step3中的选择相同；
[0137] 其中，混合信号x1(n)、x2(n)为脉搏信号x(m)(m＝1,2,3,…,M)分为长度相同的两段信号x1(i)和x2(i) 补零得到的新的信号x1(n)和x2(n)(n＝1,2,3,…,N)，N为最接近原信号x(m)长度M的2的整数次方，且N满足N≥M。
[0138] 信号n(n)(n＝1,2,3,…,N)的离散功率谱 如图13所示；
[0139] Step5、用补零后的混合信号的功率谱减去噪声信号模型的功率谱，得到有效信号的功率谱，结合去噪前混合信号的相位谱，通过变换得到时域脉搏有效信号。
[0140] 所述步骤Step5的具体步骤如下：
[0141] Step5.1、用信号x1(n)的功率谱估计 减去信号n(n)的功率谱估计得到信号s1(n)的功率谱估计 即 如图14
所示；将 记为
[0142] Step5.2、对 作逆变换，得 将 的前 项与后 项调换顺序，使后 项在前，前 项在后，得到新的序列 即
和真实自相关函数 有如下关系
[0143]
[0144] 求得 为 中间的2p+1项，即记 为
[0145] Step5.3、对 补零得到 使 的长度与 相同，为N1，且使 的前p+1项等于的前项，即 使 的后p+1项等于 的后p+1项，即
使 的中间N1-2(p+1 )项的值均等于0，即
取 然后对 进行傅里叶变换，则变换结果的实部为有
效信号s1(n)的傅里叶变换幅值的平方|S1(k)2，对变换结果的实部开平方，得到有效信号s1(n)傅里叶变换的幅值|S1(k)|； 如图15所示；
[0146] Step5.4、结合信号x1(n)的相位谱 得到信号s1(n)的傅里叶变换结果S1(k)，即 对S1(k)取傅里叶逆变换，得到有效信号s1(n)，且取信号的前I个值，即可得到包含在混合信号x1(i)中的有效信号s1(i)；如图7所示；
[0147] Step5.5、用与步骤Step5.1至Step5.4相同的方法可以得到 s2(i),将s1(i)与s2(i)首尾相接，即可得到去噪后的整段信号s(m)； s2(i)、s(m)分别如图17、图9、图5所示， 如图16所示；
[0148] 在信号x1(i)(i＝1,2,3,…,I)中截取一段信号x1_1(t)(t＝1,2,3,…,T)；其中，数据长度T＝3000，该段信号截取的是x1(i)中的前三千个采样点数据；x1_1(t)(t＝1,2,3,…,T)，如图18所示；
[0149] 在去噪后信号s1(i)(i＝1,2,3,…,I)中截取一段信号s1_1(t)(t＝1,2,3,…,T)；其中，数据长度T＝3000，该段信号截取的是s1(i)中的前三千个采样点数据，即s1_1(t)(t＝
1,2,3,…,T)就是x1_1(t)(t＝1,2,3,…,T)去噪后的结果；s1_1(t)(t＝1,2,3,…,T)，如图19所示；
[0150] 在信号x2(i)(i＝1,2,3,…,I)中截取一段信号x2_1(t)(t＝1,2,3,…,T)；其中，数据长度T＝3000，该段信号截取的是x2(i)中的前三千个采样点数据；x2_1(t)(t＝1,2,3,…,T)，如图20所示；
[0151] 在去噪后信号s2(i)(i＝1,2,3,…,I)中截取一段信号s2_1(t)(t＝1,2,3,…,T)；其中，数据长度T＝3000，该段信号截取的是s1(i)中的前三千个采样点数据，即s1_1(t)(t＝
1,2,3,…,T)就是x2_1(t)(t＝1,2,3,…,T)去噪后的结果；s2_1(t)(t＝1,2,3,…,T)，如图21所示；
[0152] 如图对比图4与图5、图6与图7、图8与图9，可以直观的看出，图5、图7、图9中所示的信号相对于图4、图6、图8中所示信号，在高度保持了原信号形状的情况下，信号更加平滑，说明很好地去除了随机噪声信号；在图18、图20所示信号是分别在x1(i)、x2(i)截取的一段信号，图20、21中所示信号是对这两段信号去噪后的结果，对比图18与图19、图20与图21中所示信号，可以清晰地看到去噪后的信号极大地保持了原有信号的形状，同时去噪后的信号更加平滑，分辨率很高，信号的拐点处更加清晰。由此可以看出，文中所给方法对脉搏波中随机噪声的去噪效果很好。
[0153] 上面结合附图对本发明的具体实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。
