一种利用先验信息的指数信号去噪方法与流程

本发明涉及指数信号的去噪方法，尤其是涉及一种利用先验信息的指数信号去噪方法。

背景技术：

在实际应用中，采集数据时往往会受到噪声污染，使得微弱的有效信号无法辨别。对于一般信号，通常采用重复采集的方式得到更高信噪比的信号，但这使采样时间延长，特别在一些大批量采集中会造成耗时过长；而对于一些特殊信号，可以利用信号特性采用一些效果更优的去噪方法。对于指数信号，比如核磁共振波谱、雷达载波信号等，通常采集得到的时域(相对频域)信号可以表示成一系列指数函数的叠加，这种特殊的信号可以建模成在频域(相对时域)上若干谱峰的线性叠加。这类信号经过构建汉克尔矩阵的线性算子作用后具有低秩特性(Y.Y.Lin and L.P.Hwang,NMR signal enhancement based on matrix property mappings[J],Journal of Magnetic Resonance,Series A,1993,103(1):109-114；X.Qu,M.Mayzel,J.-F.Cai,Z.Chen,and V.Y.Orekhov,Accelerated NMR spectroscopy with low-rank reconstruction[J],Angewandte Chemie International Edition,2015,54(3):852-854.)。基于这种特性的一种典型去噪方法是截断奇异值分解(G.Golub and W.Kahan,Calculating the singular values and pseudoinverse of a matrix[J],J.Soc.Ind.Appl.Math.B,1965,2(2):205–224)，虽然这种方法的去噪效果优于重复采集的方式，但去噪效果有限。

技术实现要素：

本发明的目的在于提供效果优良、易于操作的一种利用先验信息的指数信号去噪方法。

本发明包括以下步骤：

1)对参考指数信号构建汉克尔矩阵；

在步骤1)中，所述对参考指数信号构建汉克尔矩阵的具体方法可为：将参考指数信号记做f＝[f(1),f(2),…,f(N)]，长度记做N,其中b_j和g_j为复数，J为正整数，n为指数的阶次；通过线性算子R，将向量f构建成汉克尔矩阵F：

式(1)中的Q,P是R的两个参数，分别决定了矩阵F的行数和列数。

2)获取所需的先验信息：对步骤1)中构建的汉克尔矩阵F进行奇异值分解，得到左奇异向量矩阵U_f和奇异值s(i),i＝1,2,…,Q：

上式中V_f为右奇异向量矩阵，符号“H”表示求矩阵的共轭转置。

3)对目标指数信号构建汉克尔矩阵；

在步骤3)中，所述对目标指数信号构建汉克尔矩阵的具体方法可为：将长度同为N的指数信号记做x＝[x(1),x(2),…,x(N)]，其中c_k和z_k为复数，K为正整数，n为指数的阶次；通过线性算子R，将向量x构建成汉克尔矩阵X：

4)获得目标信号汉克尔矩阵X的奇异值：利用步骤2)中的先验信息U_f矩阵，对矩阵X进行分解，得到λ(i),i＝1,2,…,Q：

式(4)中，其中diag(·)_i表示取矩阵第i行第i列的值；其中Λ是一个对角矩阵，其对角线上第i个值为1/λ(i)。

5)设置加权的阈值：利用步骤2)中参考信号的奇异值s(i)，设置加权阈值为：其中e＞0是设置的参数，ε＝10^-6；

6)获取新的奇异值：令α(i)＝max(λ(i)-w(i),0),i＝1,2,…,Q，其中max(l₁,l₂)表示取l₁和l₂中最大数值；

7)由新奇异值得到去噪后的指数信号：利用步骤6)中的奇异值α(i)得到重组矩阵再经过逆汉克尔算子得到去噪后的指数信号

其中，符号“-1”表示求矩阵的逆，符号“T”表示求矩阵的转置；

8)数据后处理：对去噪后的指数信号进行傅立叶变换得到频谱。

本发明基于指数信号经过线性算子作用后的汉克尔矩阵的低秩特性，利用先验信息进行矩阵分解和奇异值更新，达到对信号去噪的目的。这种方法适用于在大批量的对比实验中(比如采集多种样品的核磁共振波谱图)大量节省数据采集时间。

附图说明

图1是一维核磁共振频谱的参考信号。

图2是无噪的一维核磁共振频谱信号。

图3是加噪后的一维核磁共振频谱信号。

图4是本发明方法去噪后得到的一维核磁共振频谱信号。

具体实施方式

以下具体实施例对一维核磁共振波谱进行去噪，具体步骤如下：

1)对参考指数信号构建汉克尔矩阵：图1是参考信号的频谱。将一维核磁共振波谱时域(相对于频域)参考信号记做f＝[f(1),f(2),…,f(1023)]，长度为N＝1023。通过线性算子R，将向量f构建成汉克尔矩阵F：

上式中R的两个参数Q＝512,P＝512。

2)获取所需的先验信息：对1)中构建的汉克尔矩阵F进行奇异值分解，得到左奇异向量矩阵U_f和奇异值s(i),i＝1,2,…,512：

上式中V_f为右奇异向量矩阵，符号“H”表示求矩阵的共轭转置。

3)对目标指数信号构建汉克尔矩阵：图3是对图2中无噪信号添加噪声后得到的带噪目标信号的频谱。将长度同为N的一维核磁共振波谱时域信号(相对于频域)记为x＝[x(1),x(2),…,x(1023)]。通过线性算子R，将向量x构建成汉克尔矩阵：

4)获得目标信号汉克尔矩阵X的奇异值：利用步骤2)中的先验信息U_f矩阵，对矩阵X进行分解，得到λ(i),i＝1,2,…,512：

式(9)中其中diag(·)_i表示取矩阵第i行第i列的值；其中Λ是一个对角矩阵，其对角线上第i个值为1/λ(i)。

5)设置加权的阈值：利用步骤2)中参考信号的奇异值s(i)，设置加权阈值为：其中e＝1.3，ε＝10^-6。

6)获取新的奇异值：令α(i)＝max(λ(i)-w(i),0),i＝1,2,…,512，其中max(l₁,l₂)表示取l₁和l₂中最大数值。

7)由新奇异值得到去噪后的指数信号：利用步骤6)中的奇异值α(i)得到重组矩阵再经过逆汉克尔算子得到去噪后的指数信号

其中符号“-1”表示求矩阵的逆，符号“T”表示求矩阵的转置。

8)数据后处理：对去噪后的时域信号进行傅立叶变换得到谱图(如图4所示)。