CN107729845A - 一种基于子空间特征值分解的实测频响函数降噪方法 - Google Patents

一种基于子空间特征值分解的实测频响函数降噪方法 Download PDF

Info

Publication number
CN107729845A
CN107729845A CN201710981385.8A CN201710981385A CN107729845A CN 107729845 A CN107729845 A CN 107729845A CN 201710981385 A CN201710981385 A CN 201710981385A CN 107729845 A CN107729845 A CN 107729845A
Authority
CN
China
Prior art keywords
mrow
msub
msup
msubsup
mtd
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
CN201710981385.8A
Other languages
English (en)
Other versions
CN107729845B (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.)
New Energy Automobile Group Co Ltd
Original Assignee
New Energy Automobile Group Co Ltd
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 New Energy Automobile Group Co Ltd filed Critical New Energy Automobile Group Co Ltd
Priority to CN201710981385.8A priority Critical patent/CN107729845B/zh
Publication of CN107729845A publication Critical patent/CN107729845A/zh
Application granted granted Critical
Publication of CN107729845B publication Critical patent/CN107729845B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/02—Preprocessing
    • G06F2218/04—Denoising
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G06F17/141—Discrete Fourier transforms

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Discrete Mathematics (AREA)
  • Algebra (AREA)
  • Signal Processing (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Noise Elimination (AREA)
  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)

Abstract

本发明公开了一种基于子空间特征值分解的实测频响函数降噪方法，主要包括1)初步推导最优估计器，并进行特征值分解；2)求解最优估计器，并得到降噪结果。本发明方法不同于传统的降噪方法，是在完成频响函数测试以后进一步对频响函数进行降噪。为了验证方法的有效性，设计六阶系统，并人为引入了噪声，分别利用两种算法对带噪频响函数进行了降噪分析，证明了子空间降噪算法在频响函数降噪中的有效性，最后将降噪算法应用在实测频响函数的降噪中，取得了较好的降噪效果。

Description

一种基于子空间特征值分解的实测频响函数降噪方法
技术领域：
本发明涉及一种基于子空间特征值分解的实测频响函数降噪方法，其属于频响函数降噪处理领域。
背景技术：
模态识别、有限元模型修正、故障诊断、结构损伤定位、振动噪声传递路径分析，以及动态载荷识别等分析过程都直接需要对频响函数或者频响函数矩阵进行微分或求逆运算。在子结构频响函数综合中，也需要对界面测点的频响函数矩阵进行求逆运算。求逆或者微分过程对信号噪声十分敏感。因此如何降低实测频响函数中噪声水平成为上述分析能否成功的一个关键。针对实际测量中频响函数的激励与响应存在的噪声，工程上一般采用H1、H2、H3、H4以及Hv等经典的方法，或者Hc估计等其它一些方法对实测频响函数进行估计或者说降噪处理。一般来说，工程应用中使用Hv估计方法较多。以上这些经典的方法属于在测量频响函数时为了排除噪声的影响而采用的估计方法，其中大都采用平均技术。
Sanliturk等利用奇异值分解(SVD)方法对测试中获得频响函数进行了进一步降噪处理。其基本思路是对带噪频响函数进行相空间重构，利用奇异值分解方法将重构的相空间矩阵分解为信号子空间和噪声子空间，认为前若干个较大的奇异值对应的就是信号子空间，而其它奇异值对应的是噪声子空间，并将其置为零。再利用奇异值分解的逆变换对纯净信号矩阵进行估计，最后由相空间重构的逆过程，通过平均的方法获得降噪后的频响函数。孙鑫晖等人利用同样的方法对频响函数进行了降噪。利用奇异值分解降噪是一种有效的降噪分析方法，在对信号增强与信号特征提取中有较多应用。基于奇异值分解方法实际上属于子空间降噪算法中的一种。子空间方法不同于其它时频降噪方法，如维纳滤波或者卡尔曼滤波等等。子空间降噪算法在语音增强、障诊断以及图像降噪方面均有应用。Ephraim等通过对协方差矩阵的特征值分解(EVD)，带噪语音信号进行子空间分解，并与谱减法进行了对比，证明其在降噪方面的优势。Hu等将纯净的语音信号与噪声信号的协方差同时对角化，将Ephraim方法推广到了有色噪声的降噪处理中。Hermus等人对目前基于子空间语音增强的方法进行了总结，并将各种子结构语音增强算法进行了对比，指出与谱减法相比，子空间方法减少了音乐噪声的残留，提高了语音识别的准确性。钱征文等利用子空间中的SVD分解对信号进行了降噪处理，利用FFT变换结果确定重构相空间矩阵的秩。吕勇等利用子空间降噪算法对齿轮故障信号进行了分析与降噪，提取出相应的故障特征信号。
综上，H1-H4以及Hv等经典的方法，均是在完成激励与响应点响应测量之后，在获得频响函数之前通过互谱平均的方法对频响函数进行估计，以获得信噪比较高的频响函数。但通常利用这些经典的方法估计后获得的频响函数依然会有噪声残留，这些残留噪声依然会影响频响函数逆运算或微分运算的精度。
发明内容：
本发明是为了解决上述现有技术存在的问题而提供一种基于子空间特征值分解的实测频响函数降噪方法，能够有效地除去频响函数的背景噪声。
本发明采用如下技术方案：一种基于子空间特征值分解的实测频响函数降噪方法，包括如下步骤：
步骤1、频响函数的离散傅里叶反变换，子空间算法在时域内进行，根据傅里叶反变换将频响函数变为时域内的脉冲响应函数对于带噪的频响函数，h为脉冲响应函数，H与h均为离散观测序列，n为脉冲响应序列的长度,据离散傅里叶反变换
将离散傅里叶反变换后的带噪脉冲响应函数写成由纯净的脉冲响应函数h'k加上噪声nk构成
hk＝h'k+nk,0≤k≤n-1 (2)
假定频域噪声信号Ns的实部与虚部独立，服从均值为零的高斯分布，且方差相等即：
Real(Ns)～G(0,σ2),Imag(Ns)～G(0,σ2) (3)
其中Ns是nk的复数域表达，根据离散傅里叶变换对(DFT)：
将nk展开成实部与虚部
利用欧拉公式展开，经推导可证明nk也服从高斯分布，期望与方差分别为:
E[nk]＝0 (6)
即对于频域带噪信号，若其噪声信号的实部与虚部独立，并服从均值为零的高斯分布，则经过傅里叶反变换后的脉冲响应函数同样服从高斯分布，期望与方差如上式；
步骤2、构造线性滤波器L，将由纯净的脉冲响应函数h'k加上噪声nk构成的带噪脉冲响应函数成向量的形式：
h＝h′+n (8)
上式中向量的长度均为n×1，假定h′的估计可通过h的线性滤波器L构成：
将残余误差写成由频响函数失真误差εh′和残余噪声误差εn两个部分构成
其中，频响函数失真误差定义为：
残余噪声误差定义为：
由上述式(6)、(7)、以及aTa＝trace(aaT)推导出频响函数失真能量为：
残余噪声能量为：
其中，R*为协方差矩阵，trace为矩阵的迹，由上述推导可得到以下时域约束最优问题用于求解最优线性估计器：
其中,δ为正常数，即满足残余噪声能量小于δ时，频响函数失真最小，引入拉格朗日乘子μ，μ≥0，构造拉格朗日乘式：
根据KKT条件:最终可得到最优估计器：
步骤3、特征值分解，假定嵌入维数为m，延迟时间为τ，利用延迟坐标法重构脉冲响应函数的相空间矩阵，重构的相空间矩阵Y写成：
对降噪分析而言，直接取延迟时间τ＝1，重构的相空间矩阵Y写成：
其中，Y为Toeplitz矩阵，假设频响函数中存在加性噪声，则相空间矩阵Y可写成：
Y＝X+N (20)
其中，Y、X和N分别表示原始带噪、纯净信号和噪声信号的脉冲响应函数的重构矩阵，带噪脉冲响应函数的协方差可写成：
Rh＝Rh′+Rn (21)
采用下式对Rh进行估计：
步骤4、用Rh对Rh′进行估计，Rh为m×m矩阵,将纯净脉冲响应函数的协方差矩阵Rh′进行EVD分解：
将式(23)带入式(21)：
假定Rh′的秩为r，其中r＜m，则Λh′中的前r个对角元素不为零，即
将式(25)带入(24)，则原始脉冲响应的协方差矩阵Rh可写成分块形式：
直接对通过估计获得的Rh进行EVD分解，并按(26)进行分块：
对比式(26)、(27)，可知：
将式(28)带入式(25)，可知：
步骤5、得出最优估计器，根据式(17)、(23)以及(29)，最终得最优估计器LTDC
当μ＝1，式(30)即为所谓线性最小二乘误差(LMMSE)估计器：
当μ＝0，即为所谓最小二乘(LS)估计器：
步骤6、求解纯净脉冲响应函数，根据式(30)可知纯净脉冲响应函数估计表示成
步骤7、求解降噪后频响函数，在获得纯净脉冲响应函数估计后，再利用傅里叶变换即可得到降噪后频响函数；
步骤8、求解降噪结果，依据特征值分析结果，对比保留不同特征值个数对应的降噪效果，确定协方差特征值个数，估算其后验信噪比SNRdB，并取拉格朗日乘子μ，依据最优估计器LTDC求解得到降噪结果。
本发明具有如下有益效果：考虑实测频响函数噪声问题，本发明提出了基于子空间特征值分解的实测频响函数降噪方法，该方法不同于传统的降噪方法，是在完成频响函数测试以后进一步对频响函数进行降噪。为了验证方法的有效性，设计六阶系统，并人为引入了噪声，分别利用两种算法对带噪频响函数进行了降噪分析，证明了子空间降噪算法在频响函数降噪中的有效性，最后将降噪算法应用在实测频响函数的降噪中，取得了较好的降噪效果。
附图说明：
图1为纯净与带噪脉冲响应函数。
图2为纯净与带噪频响函数幅值与相位对比。
图3为协方差矩阵特征值。
图4为不同μ的降噪效果幅值与相位对比。
图5为SVD与EVD降噪效果对比。
图6为实测频响函数Hs2x,8z幅值与相位。
图7为脉冲响应函数。
图8为特征值分析结果。
图9为降噪前后频响函数Hs2x,8z幅值与相位对比。
具体实施方式：
下面结合附图对本发明作进一步的说明。
本发明基于子空间特征值分解的实测频响函数降噪方法，主要包括：
1)初步推导最优估计器，并进行特征值分解；
2)求解最优估计器，并得到降噪结果。
进一步地，具体包括如下步骤：
步骤1：频响函数的离散傅里叶反变换。子空间算法一般在时域内进行，根据傅里叶反变换将频响函数变为时域内的脉冲响应函数对于带噪的频响函数，h为脉冲响应函数，这里H与h均为离散观测序列，n为脉冲响应序列的长度,据离散傅里叶反变换
将离散傅里叶反变换后的带噪脉冲响应函数写成由纯净的脉冲响应函数h'k加上噪声nk构成
hk＝h'k+nk,0≤k≤n-1 (2)
假定频域噪声信号Ns的实部与虚部独立，服从均值为零的高斯分布，且方差相等即：
Real(Ns)～G(0,σ2),Imag(Ns)～G(0,σ2) (3)
其中Ns是nk的复数域表达，根据离散傅里叶变换对(DFT)：
将nk展开成实部与虚部
利用欧拉公式展开，经推导可证明nk也服从高斯分布，期望与方差分别为:
E[nk]＝0 (6)
即对于频域带噪信号，若其噪声信号的实部与虚部独立，并服从均值为零的高斯分布，则经过傅里叶反变换后的脉冲响应函数同样服从高斯分布，期望与方差如上式。
步骤2：构造线性滤波器L。将由纯净的脉冲响应函数h'k加上噪声nk构成的带噪脉冲响应函数成向量的形式：
h＝h′+n (8)
上式中向量的长度均为n×1。假定h′的估计可通过h的线性滤波器L构成：
将残余误差可写成由频响函数失真误差εh′和残余噪声误差εn两个部分构成
其目的是保证残余误差低于一定阈值的同时使频响函数失真度最小。显然该方式是频响函数失真和残余误差程度之间的一种折中。其中，
频响函数失真误差定义为：
残余噪声误差定义为：
由上述式(6)、(7)、以及aTa＝trace(aaT)可推导出频响函数失真能量为：
残余噪声能量为：
其中，R*为协方差矩阵，trace为矩阵的迹。
由上述推导可得到以下时域约束最优问题用于求解最优线性估计器：
其中,δ为正常数，即满足残余噪声能量小于δ时，频响函数失真最小。引入拉格朗日乘子μ，μ≥0，构造拉格朗日乘式：
根据KKT条件:最终可得到最优估计器：
可见，如果已知Rh′，即可得到最优估计器LTDC，但通常来说Rh′是未知的，需要利用Rh进行估计。
步骤3：特征值分解。假定嵌入维数为m，延迟时间为τ，利用延迟坐标法重构脉冲响应函数的相空间矩阵。重构的相空间矩阵Y可以写成：
对降噪分析而言，一般直接取延迟时间τ＝1，重构的相空间矩阵Y可以写成：
其中，Y为Toeplitz矩阵，当然也可以将脉冲响应序列重构成Hankel矩阵形式。假设频响函数中存在加性噪声，则相空间矩阵Y可写成：
Y＝X+N (20)
其中，Y、X和N分别表示原始带噪、纯净信号和噪声信号的脉冲响应函数的重构矩阵。带噪脉冲响应函数的协方差可写成：
Rh＝Rh′+Rn (21)
可采用下式对Rh进行估计：
步骤4：用Rh对Rh′进行估计。Rh为m×m矩阵,将纯净脉冲响应函数的协方差矩阵Rh′进行EVD分解：
将式(23)带入式(21)：
假定Rh′的秩为r，其中r＜m，则Λh′中的前r个对角元素不为零，即
将式(25)带入(24)，则原始脉冲响应的协方差矩阵Rh可写成分块形式：
直接对通过估计获得的Rh进行EVD分解，并按(26)进行分块：
对比式(26)、(27)，可知：
将式(28)带入式(25)，可知：
步骤5：得出最优估计器。根据式(17)、(23)以及(29)，最终可得最优估计器LTDC
事实上，当μ＝1，式(30)即为所谓线性最小二乘误差(LMMSE)估计器：
当μ＝0，即为所谓最小二乘(LS)估计器：
步骤6：求解纯净脉冲响应函数。根据式(30)可知纯净脉冲响应函数估计可表示成
步骤7：求解降噪后频响函数。在获得纯净脉冲响应函数估计后，再利用傅里叶变换即可得到降噪后频响函数。
步骤8：求解降噪结果。依据特征值分析结果，对比保留不同特征值个数对应的降噪效果，确定协方差特征值个数，估算其后验信噪比SNRdB，并取拉格朗日乘子μ，依据最优估计器LTDC求解得到降噪结果。
算例一：
为验证基于子空间特征值分解的算法在频响函数降噪中的效果，设计一个六阶系统，其传递函数为其中s为Laplace算子。系统参数设定如下:m1＝1，m2＝0.5，m3＝0.8，k1＝40π2，k2＝400π2，k3＝4000π2， 给定采样率为100Hz，数据样本长度为1000，令传递函数中s＝jω，计算纯净频响函数的单位脉冲响应，并叠加“加性”高斯白噪声，Matlab的命令为awgn(x,10,'measured')，信噪比SNR为10dB。
加入噪声前、后的脉冲响应函数如图1。
经过傅里叶变换后的频响函数幅值与相位如图2。
步骤1：特征值分解。利用式(22)估计脉冲响应函数重构相空间矩阵的协方差矩阵，并进行特征值分解。图3给出了前50个特征值的分析结果。
步骤2：最优估计器与降噪结果。根据图3，取协方差特征值个数为6。注意式(30)，最优估计器LTDC与拉格朗日乘子μ有关。图4给出μ＝0～10时降噪结果。
由图4可知，显然降噪后频响函数无论是幅值与相位，其精度大幅度提高，特别是对频响函数的反共振点的估计。但降噪效果与μ的大小并非是单调的，观察发现当μ∈(0,5)时，降噪效果最佳。目前关于μ取值尚无统一的计算公式，这里给出语音增强算法中μ取值的一种估计方法供参考。
其中,μ＝4.2，c＝6.25，SNRdB为后验信噪比估计。这里SNRdB约为10dB，计算可得μ＝2.6，与前述分析的μ∈(0,5)的最优取值范围是一致的。
步骤3：降噪效果对比。为了检验在低信噪比情况下基于EVD与SVD两种算法的降噪效果，将信噪比下调至5dB。
图5为采用SVD_LS方法与EVD(μ＝3)两种方法降噪后的结果。其中，SVD_LS为基于子空间SVD最小二乘意义下的降噪算法。
显然即使在信噪比达到5dB时两种方法依然具有较好的降噪效果，而从降噪后的相位对比来看，EVD方法降噪效果比SVD_LS方法更优。
算例二：
以某微车的实测跨点频响函数为例，利用EVD方法的时域约束估计方法对实测频响函数进行降噪。其中频响函数通过锤击实验法得到，采样率为1024Hz，频率分辨率为0.125Hz，频响函数的估计方法采用H1估计。
图6给出了车身左前悬架+z激励下驾驶室内右前导轨+x测点响应的实测频响函数。
步骤1：特征值分解。图7给出了经过傅里叶反变换得到的实测频响函数的脉冲响应函数。构造Topelitz的频响函数相空间矩阵Y，并对Y进行特征值分解.
步骤2：最优估计器与降噪结果及对比。图8给出了脉冲响应函数重构相空间矩阵的协方差矩阵的前100个特征值的分析结果。根据图8特征值分析结果，并对比保留不同特征值个数对应的降噪效果，最终确定协方差特征值个数为30，经过估算其后验信噪比SNRdB约为16dB，最终取拉格朗日乘子μ＝1.6，并由式(30)，确定最优估计器LTDC。
图9给出了TDC降噪结果，由图可知采用EVD方法的时域约束估计方法降噪结果。显然，结果降噪后能明显提高频响函数的精度，特别是对反共振点的改善。
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下还可以作出若干改进，这些改进也应视为本发明的保护范围。

Claims (1)

1.一种基于子空间特征值分解的实测频响函数降噪方法，其特征在于：包括如下步骤：
步骤1、频响函数的离散傅里叶反变换，子空间算法在时域内进行，根据傅里叶反变换将频响函数变为时域内的脉冲响应函数对于带噪的频响函数，h为脉冲响应函数，H与h均为离散观测序列，n为脉冲响应序列的长度,据离散傅里叶反变换
<mrow> <msub> <mi>h</mi> <mi>k</mi> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mi>n</mi> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>s</mi> <mo>=</mo> <mn>0</mn> </mrow> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </munderover> <msub> <mi>H</mi> <mi>s</mi> </msub> <msup> <mi>e</mi> <mrow> <mi>j</mi> <mfrac> <mrow> <mn>2</mn> <mi>&amp;pi;</mi> <mi>s</mi> <mi>k</mi> </mrow> <mi>n</mi> </mfrac> </mrow> </msup> <mo>,</mo> <mn>0</mn> <mo>&amp;le;</mo> <mi>k</mi> <mo>&amp;le;</mo> <mi>n</mi> <mo>-</mo> <mn>1</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>
将离散傅里叶反变换后的带噪脉冲响应函数写成由纯净的脉冲响应函数h'k加上噪声nk构成
hk＝h'k+nk,0≤k≤n-1 (2)
假定频域噪声信号Ns的实部与虚部独立，服从均值为零的高斯分布，且方差相等即：
Real(Ns)～G(0,σ2),Imag(Ns)～G(0,σ2) (3)
其中Ns是nk的复数域表达，根据离散傅里叶变换对(DFT)：
<mrow> <msub> <mi>N</mi> <mi>s</mi> </msub> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>0</mn> </mrow> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </munderover> <msub> <mi>n</mi> <mi>k</mi> </msub> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>j</mi> <mn>2</mn> <mi>&amp;pi;</mi> <mi>k</mi> <mi>s</mi> <mo>/</mo> <mi>n</mi> </mrow> </msup> <mo>,</mo> <msub> <mi>n</mi> <mi>k</mi> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mi>n</mi> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>0</mn> </mrow> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </munderover> <msub> <mi>N</mi> <mi>s</mi> </msub> <msup> <mi>e</mi> <mrow> <mo>+</mo> <mi>j</mi> <mn>2</mn> <mi>&amp;pi;</mi> <mi>k</mi> <mi>s</mi> <mo>/</mo> <mi>n</mi> </mrow> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow>
将nk展开成实部与虚部
<mrow> <msub> <mi>n</mi> <mi>k</mi> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mi>n</mi> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>0</mn> </mrow> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </munderover> <msub> <mi>N</mi> <mi>s</mi> </msub> <msup> <mi>e</mi> <mrow> <mo>+</mo> <mi>j</mi> <mn>2</mn> <mi>&amp;pi;</mi> <mi>k</mi> <mi>s</mi> <mo>/</mo> <mi>n</mi> </mrow> </msup> <mo>=</mo> <mfrac> <mn>1</mn> <mi>n</mi> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>0</mn> </mrow> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </munderover> <mi>Re</mi> <mi>a</mi> <mi>l</mi> <mrow> <mo>(</mo> <msub> <mi>N</mi> <mi>s</mi> </msub> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mrow> <mo>+</mo> <mi>j</mi> <mn>2</mn> <mi>&amp;pi;</mi> <mi>k</mi> <mi>s</mi> <mo>/</mo> <mi>n</mi> </mrow> </msup> <mo>+</mo> <mfrac> <mn>1</mn> <mi>n</mi> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>0</mn> </mrow> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </munderover> <mi>Im</mi> <mi>a</mi> <mi>g</mi> <mrow> <mo>(</mo> <msub> <mi>N</mi> <mi>s</mi> </msub> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mrow> <mo>+</mo> <mi>j</mi> <mn>2</mn> <mi>&amp;pi;</mi> <mi>k</mi> <mi>s</mi> <mo>/</mo> <mi>n</mi> </mrow> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow>
利用欧拉公式展开，经推导可证明nk也服从高斯分布，期望与方差分别为:
E[nk]＝0 (6)
<mrow> <mi>V</mi> <mi>a</mi> <mi>r</mi> <mo>&amp;lsqb;</mo> <msub> <mi>n</mi> <mi>k</mi> </msub> <mo>&amp;rsqb;</mo> <mo>=</mo> <mfrac> <mn>1</mn> <msup> <mi>n</mi> <mn>2</mn> </msup> </mfrac> <msup> <mi>n&amp;sigma;</mi> <mn>2</mn> </msup> <mo>+</mo> <mfrac> <mn>1</mn> <msup> <mi>n</mi> <mn>2</mn> </msup> </mfrac> <msup> <mi>n&amp;sigma;</mi> <mn>2</mn> </msup> <mo>=</mo> <mfrac> <mn>2</mn> <mi>n</mi> </mfrac> <msup> <mi>&amp;sigma;</mi> <mn>2</mn> </msup> <mo>=</mo> <msubsup> <mi>&amp;sigma;</mi> <mrow> <mi>n</mi> <mi>o</mi> <mi>i</mi> <mi>s</mi> <mi>e</mi> </mrow> <mn>2</mn> </msubsup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow>
即对于频域带噪信号，若其噪声信号的实部与虚部独立，并服从均值为零的高斯分布，则经过傅里叶反变换后的脉冲响应函数同样服从高斯分布，期望与方差如上式；
步骤2、构造线性滤波器L，将由纯净的脉冲响应函数h'k加上噪声nk构成的带噪脉冲响应函数成向量的形式：
h＝h′+n (8)
上式中向量的长度均为n×1，假定h′的估计可通过h的线性滤波器L构成：
<mrow> <mover> <mi>h</mi> <mo>^</mo> </mover> <mo>=</mo> <mi>L</mi> <mi>h</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> </mrow>
将残余误差写成由频响函数失真误差εh′和残余噪声误差εn两个部分构成
其中，频响函数失真误差定义为：
<mrow> <msub> <mi>&amp;epsiv;</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> </msub> <mover> <mo>=</mo> <mi>&amp;Delta;</mi> </mover> <mrow> <mo>(</mo> <mi>L</mi> <mo>-</mo> <mi>E</mi> <mo>)</mo> </mrow> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>11</mn> <mo>)</mo> </mrow> </mrow>
残余噪声误差定义为：
<mrow> <msub> <mi>&amp;epsiv;</mi> <mi>n</mi> </msub> <mover> <mo>=</mo> <mi>&amp;Delta;</mi> </mover> <mi>L</mi> <mi>n</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>12</mn> <mo>)</mo> </mrow> </mrow>
由上述式(6)、(7)、以及aTa＝trace(aaT)推导出频响函数失真能量为：
<mrow> <mover> <msubsup> <mi>&amp;epsiv;</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mn>2</mn> </msubsup> <mo>&amp;OverBar;</mo> </mover> <mo>=</mo> <mi>E</mi> <mrow> <mo>(</mo> <msubsup> <mi>&amp;epsiv;</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mi>T</mi> </msubsup> <msub> <mi>&amp;epsiv;</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mi>t</mi> <mi>r</mi> <mi>a</mi> <mi>c</mi> <mi>e</mi> <mrow> <mo>(</mo> <mi>E</mi> <mo>&amp;lsqb;</mo> <msub> <mi>&amp;epsiv;</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> </msub> <msubsup> <mi>&amp;epsiv;</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mi>T</mi> </msubsup> <mo>&amp;rsqb;</mo> <mo>)</mo> </mrow> <mo>=</mo> <mi>t</mi> <mi>r</mi> <mi>a</mi> <mi>c</mi> <mi>e</mi> <mo>&amp;lsqb;</mo> <mrow> <mo>(</mo> <mi>L</mi> <mo>-</mo> <mi>E</mi> <mo>)</mo> </mrow> <msub> <mi>R</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> </msub> <msup> <mrow> <mo>(</mo> <mi>L</mi> <mo>-</mo> <mi>E</mi> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mo>&amp;rsqb;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>13</mn> <mo>)</mo> </mrow> </mrow>
残余噪声能量为：
<mrow> <mover> <msubsup> <mi>&amp;epsiv;</mi> <mi>n</mi> <mn>2</mn> </msubsup> <mo>&amp;OverBar;</mo> </mover> <mo>=</mo> <mi>E</mi> <mrow> <mo>(</mo> <msubsup> <mi>&amp;epsiv;</mi> <mi>n</mi> <mi>T</mi> </msubsup> <msub> <mi>&amp;epsiv;</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mi>t</mi> <mi>r</mi> <mi>a</mi> <mi>c</mi> <mi>e</mi> <mrow> <mo>(</mo> <mi>E</mi> <mo>&amp;lsqb;</mo> <msub> <mi>&amp;epsiv;</mi> <mi>n</mi> </msub> <msubsup> <mi>&amp;epsiv;</mi> <mi>n</mi> <mi>T</mi> </msubsup> <mo>&amp;rsqb;</mo> <mo>)</mo> </mrow> <mo>=</mo> <mi>t</mi> <mi>r</mi> <mi>a</mi> <mi>c</mi> <mi>e</mi> <mrow> <mo>(</mo> <msub> <mi>LR</mi> <mi>n</mi> </msub> <msup> <mi>L</mi> <mi>T</mi> </msup> <mo>)</mo> </mrow> <mo>=</mo> <msubsup> <mi>&amp;sigma;</mi> <mrow> <mi>n</mi> <mi>o</mi> <mi>i</mi> <mi>s</mi> <mi>e</mi> </mrow> <mn>2</mn> </msubsup> <mi>t</mi> <mi>r</mi> <mi>a</mi> <mi>c</mi> <mi>e</mi> <mrow> <mo>(</mo> <msup> <mi>LL</mi> <mi>T</mi> </msup> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>14</mn> <mo>)</mo> </mrow> </mrow>
其中，R*为协方差矩阵，trace为矩阵的迹，由上述推导可得到以下时域约束最优问题用于求解最优线性估计器：
<mrow> <munder> <mrow> <mi>m</mi> <mi>i</mi> <mi>n</mi> </mrow> <mi>L</mi> </munder> <mover> <msubsup> <mi>&amp;epsiv;</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mn>2</mn> </msubsup> <mo>&amp;OverBar;</mo> </mover> </mrow>
<mrow> <mi>s</mi> <mo>.</mo> <mi>t</mi> <mo>.</mo> <mfrac> <mn>1</mn> <mi>n</mi> </mfrac> <mover> <msubsup> <mi>&amp;epsiv;</mi> <mi>n</mi> <mn>2</mn> </msubsup> <mo>&amp;OverBar;</mo> </mover> <mo>&amp;le;</mo> <mi>&amp;delta;</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>15</mn> <mo>)</mo> </mrow> </mrow>
其中,δ为正常数，即满足残余噪声能量小于δ时，频响函数失真最小，引入拉格朗日乘子μ，μ≥0，构造拉格朗日乘式：
根据KKT条件:以及最终可得到最优估计器：
<mrow> <msub> <mi>L</mi> <mrow> <mi>T</mi> <mi>D</mi> <mi>C</mi> </mrow> </msub> <mo>=</mo> <msub> <mi>R</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> </msub> <msup> <mrow> <mo>(</mo> <msub> <mi>R</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> </msub> <mo>+</mo> <msub> <mi>&amp;mu;R</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mo>=</mo> <msub> <mi>R</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> </msub> <msup> <mrow> <mo>(</mo> <msub> <mi>R</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> </msub> <mo>+</mo> <msubsup> <mi>&amp;mu;&amp;sigma;</mi> <mrow> <mi>n</mi> <mi>o</mi> <mi>i</mi> <mi>s</mi> <mi>e</mi> </mrow> <mn>2</mn> </msubsup> <msub> <mi>E</mi> <mi>m</mi> </msub> <mo>)</mo> </mrow> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>17</mn> <mo>)</mo> </mrow> </mrow>
步骤3、特征值分解，假定嵌入维数为m，延迟时间为τ，利用延迟坐标法重构脉冲响应函数的相空间矩阵，重构的相空间矩阵Y写成：
对降噪分析而言，直接取延迟时间τ＝1，重构的相空间矩阵Y写成：
其中，Y为Toeplitz矩阵，假设频响函数中存在加性噪声，则相空间矩阵Y可写成：
Y＝X+N (20)
其中，Y、X和N分别表示原始带噪、纯净信号和噪声信号的脉冲响应函数的重构矩阵，带噪脉冲响应函数的协方差可写成：
Rh＝Rh′+Rn (21)
采用下式对Rh进行估计：
<mrow> <msub> <mi>R</mi> <mi>h</mi> </msub> <mo>&amp;ap;</mo> <mfrac> <mn>1</mn> <mrow> <mi>n</mi> <mo>-</mo> <mi>m</mi> <mo>+</mo> <mn>1</mn> </mrow> </mfrac> <msup> <mi>Y</mi> <mi>T</mi> </msup> <mi>Y</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>22</mn> <mo>)</mo> </mrow> </mrow>
步骤4、用Rh对Rh′进行估计，Rh为m×m矩阵,将纯净脉冲响应函数的协方差矩阵Rh′进行EVD分解：
<mrow> <msub> <mi>R</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> </msub> <mo>=</mo> <msub> <mi>U</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> </msub> <msub> <mi>&amp;Lambda;</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> </msub> <msubsup> <mi>U</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mi>T</mi> </msubsup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>23</mn> <mo>)</mo> </mrow> </mrow>
将式(23)带入式(21)：
<mrow> <msub> <mi>R</mi> <mi>h</mi> </msub> <mo>=</mo> <msub> <mi>R</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> </msub> <mo>+</mo> <msub> <mi>R</mi> <mi>n</mi> </msub> <mo>=</mo> <msub> <mi>U</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> </msub> <msub> <mi>&amp;Lambda;</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> </msub> <msubsup> <mi>U</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mi>T</mi> </msubsup> <mo>+</mo> <msubsup> <mi>&amp;sigma;</mi> <mrow> <mi>n</mi> <mi>o</mi> <mi>i</mi> <mi>s</mi> <mi>e</mi> </mrow> <mn>2</mn> </msubsup> <msub> <mi>U</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> </msub> <msub> <mi>E</mi> <mi>m</mi> </msub> <msubsup> <mi>U</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mi>T</mi> </msubsup> <mo>=</mo> <msub> <mi>U</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> </msub> <mrow> <mo>(</mo> <msub> <mi>&amp;Lambda;</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> </msub> <mo>+</mo> <msubsup> <mi>&amp;sigma;</mi> <mrow> <mi>n</mi> <mi>o</mi> <mi>i</mi> <mi>s</mi> <mi>e</mi> </mrow> <mn>2</mn> </msubsup> <msub> <mi>E</mi> <mi>m</mi> </msub> <mo>)</mo> </mrow> <msubsup> <mi>U</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mi>T</mi> </msubsup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>24</mn> <mo>)</mo> </mrow> </mrow>
假定Rh′的秩为r，其中r＜m，则Λh′中的前r个对角元素不为零，即
<mrow> <msub> <mi>R</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> </msub> <mo>=</mo> <mo>&amp;lsqb;</mo> <mtable> <mtr> <mtd> <msub> <mi>U</mi> <mrow> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mn>1</mn> </mrow> </msub> </mtd> <mtd> <msub> <mi>U</mi> <mrow> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mn>2</mn> </mrow> </msub> </mtd> </mtr> </mtable> <mo>&amp;rsqb;</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>&amp;Lambda;</mi> <mrow> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mn>1</mn> </mrow> </msub> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> </mtable> </mfenced> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msubsup> <mi>U</mi> <mrow> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mn>1</mn> </mrow> <mi>T</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <msubsup> <mi>U</mi> <mrow> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mn>2</mn> </mrow> <mi>T</mi> </msubsup> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <msub> <mi>U</mi> <mrow> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mn>1</mn> </mrow> </msub> <msub> <mi>&amp;Lambda;</mi> <mrow> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mn>1</mn> </mrow> </msub> <msubsup> <mi>U</mi> <mrow> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mn>1</mn> </mrow> <mi>T</mi> </msubsup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>25</mn> <mo>)</mo> </mrow> </mrow>
将式(25)带入(24)，则原始脉冲响应的协方差矩阵Rh可写成分块形式：
<mrow> <msub> <mi>R</mi> <mi>h</mi> </msub> <mo>=</mo> <mo>&amp;lsqb;</mo> <mtable> <mtr> <mtd> <msub> <mi>U</mi> <mrow> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mn>1</mn> </mrow> </msub> </mtd> <mtd> <msub> <mi>U</mi> <mrow> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mn>2</mn> </mrow> </msub> </mtd> </mtr> </mtable> <mo>&amp;rsqb;</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <msub> <mi>A</mi> <mrow> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mn>1</mn> <mo>+</mo> </mrow> </msub> <msubsup> <mi>&amp;sigma;</mi> <mrow> <mi>n</mi> <mi>o</mi> <mi>i</mi> <mi>s</mi> <mi>e</mi> </mrow> <mn>2</mn> </msubsup> <msub> <mi>E</mi> <mi>r</mi> </msub> </mrow> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mrow> <msubsup> <mi>&amp;sigma;</mi> <mrow> <mi>n</mi> <mi>o</mi> <mi>i</mi> <mi>s</mi> <mi>e</mi> </mrow> <mn>2</mn> </msubsup> <msub> <mi>E</mi> <mrow> <mi>m</mi> <mo>-</mo> <mi>r</mi> </mrow> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msubsup> <mi>U</mi> <mrow> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mn>1</mn> </mrow> <mi>T</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <msubsup> <mi>U</mi> <mrow> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mn>2</mn> </mrow> <mi>T</mi> </msubsup> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>26</mn> <mo>)</mo> </mrow> </mrow>
直接对通过估计获得的Rh进行EVD分解，并按(26)进行分块：
<mrow> <msub> <mi>R</mi> <mi>h</mi> </msub> <mo>=</mo> <mo>&amp;lsqb;</mo> <mtable> <mtr> <mtd> <msub> <mi>U</mi> <mrow> <mi>h</mi> <mn>1</mn> </mrow> </msub> </mtd> <mtd> <msub> <mi>U</mi> <mrow> <mi>h</mi> <mn>2</mn> </mrow> </msub> </mtd> </mtr> </mtable> <mo>&amp;rsqb;</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>&amp;Lambda;</mi> <mrow> <mi>h</mi> <mn>1</mn> </mrow> </msub> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <msub> <mi>&amp;Lambda;</mi> <mrow> <mi>h</mi> <mn>2</mn> </mrow> </msub> </mtd> </mtr> </mtable> </mfenced> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msubsup> <mi>U</mi> <mrow> <mi>h</mi> <mn>1</mn> </mrow> <mi>T</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <msubsup> <mi>U</mi> <mrow> <mi>h</mi> <mn>2</mn> </mrow> <mi>T</mi> </msubsup> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <mo>&amp;lsqb;</mo> <mtable> <mtr> <mtd> <msub> <mi>U</mi> <mrow> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mn>1</mn> </mrow> </msub> </mtd> <mtd> <msub> <mi>U</mi> <mrow> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mn>2</mn> </mrow> </msub> </mtd> </mtr> </mtable> <mo>&amp;rsqb;</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>&amp;Lambda;</mi> <mrow> <mi>h</mi> <mn>1</mn> </mrow> </msub> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <msub> <mi>&amp;Lambda;</mi> <mrow> <mi>h</mi> <mn>2</mn> </mrow> </msub> </mtd> </mtr> </mtable> </mfenced> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msubsup> <mi>U</mi> <mrow> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mn>1</mn> </mrow> <mi>T</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <msubsup> <mi>U</mi> <mrow> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mn>2</mn> </mrow> <mi>T</mi> </msubsup> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>27</mn> <mo>)</mo> </mrow> </mrow>
对比式(26)、(27)，可知：
<mrow> <msub> <mi>R</mi> <mi>h</mi> </msub> <msub> <mi>&amp;Lambda;</mi> <mrow> <mi>h</mi> <mn>1</mn> </mrow> </msub> <mo>=</mo> <msub> <mi>&amp;Lambda;</mi> <mrow> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mn>1</mn> </mrow> </msub> <mo>+</mo> <msubsup> <mi>&amp;sigma;</mi> <mrow> <mi>n</mi> <mi>o</mi> <mi>i</mi> <mi>s</mi> <mi>e</mi> </mrow> <mn>2</mn> </msubsup> <msub> <mi>E</mi> <mi>r</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>28</mn> <mo>)</mo> </mrow> </mrow>
将式(28)带入式(25)，可知：
<mrow> <msub> <mi>R</mi> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> </msub> <mo>=</mo> <msub> <mi>U</mi> <mrow> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>&amp;Lambda;</mi> <mrow> <mi>h</mi> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msubsup> <mi>&amp;sigma;</mi> <mrow> <mi>n</mi> <mi>o</mi> <mi>i</mi> <mi>s</mi> <mi>e</mi> </mrow> <mn>2</mn> </msubsup> <msub> <mi>E</mi> <mi>r</mi> </msub> <mo>)</mo> </mrow> <msubsup> <mi>U</mi> <mrow> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mn>1</mn> </mrow> <mi>T</mi> </msubsup> <mo>=</mo> <msub> <mi>U</mi> <mrow> <mi>h</mi> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>&amp;Lambda;</mi> <mrow> <mi>h</mi> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msubsup> <mi>&amp;sigma;</mi> <mrow> <mi>n</mi> <mi>o</mi> <mi>i</mi> <mi>s</mi> <mi>e</mi> </mrow> <mn>2</mn> </msubsup> <msub> <mi>E</mi> <mi>r</mi> </msub> <mo>)</mo> </mrow> <msubsup> <mi>U</mi> <mrow> <mi>h</mi> <mn>1</mn> </mrow> <mi>T</mi> </msubsup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>29</mn> <mo>)</mo> </mrow> </mrow>
步骤5、得出最优估计器，根据式(17)、(23)以及(29)，最终得最优估计器LTDC
<mrow> <msub> <mi>L</mi> <mrow> <mi>T</mi> <mi>D</mi> <mi>C</mi> </mrow> </msub> <mo>=</mo> <msub> <mi>U</mi> <mrow> <mi>h</mi> <mn>1</mn> </mrow> </msub> <mfrac> <msub> <mi>&amp;Lambda;</mi> <mrow> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mn>1</mn> </mrow> </msub> <mrow> <msub> <mi>&amp;Lambda;</mi> <mrow> <msup> <mi>h</mi> <mo>&amp;prime;</mo> </msup> <mn>1</mn> </mrow> </msub> <mo>+</mo> <msubsup> <mi>&amp;mu;&amp;sigma;</mi> <mrow> <mi>n</mi> <mi>o</mi> <mi>i</mi> <mi>s</mi> <mi>e</mi> </mrow> <mn>2</mn> </msubsup> <msub> <mi>E</mi> <mi>r</mi> </msub> </mrow> </mfrac> <msubsup> <mi>U</mi> <mrow> <mi>h</mi> <mn>1</mn> </mrow> <mi>T</mi> </msubsup> <mo>=</mo> <msub> <mi>U</mi> <mrow> <mi>h</mi> <mn>1</mn> </mrow> </msub> <mfrac> <mrow> <msub> <mi>&amp;Lambda;</mi> <mrow> <mi>h</mi> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msubsup> <mi>&amp;sigma;</mi> <mrow> <mi>n</mi> <mi>o</mi> <mi>i</mi> <mi>s</mi> <mi>e</mi> </mrow> <mn>2</mn> </msubsup> <msub> <mi>E</mi> <mi>r</mi> </msub> </mrow> <mrow> <msub> <mi>&amp;Lambda;</mi> <mrow> <mi>h</mi> <mn>1</mn> </mrow> </msub> <mo>+</mo> <mrow> <mo>(</mo> <mi>&amp;mu;</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <msubsup> <mi>&amp;sigma;</mi> <mrow> <mi>n</mi> <mi>o</mi> <mi>i</mi> <mi>s</mi> <mi>e</mi> </mrow> <mn>2</mn> </msubsup> <msub> <mi>E</mi> <mi>r</mi> </msub> </mrow> </mfrac> <msubsup> <mi>U</mi> <mrow> <mi>h</mi> <mn>1</mn> </mrow> <mi>T</mi> </msubsup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>30</mn> <mo>)</mo> </mrow> </mrow>
当μ＝1，式(30)即为所谓线性最小二乘误差(LMMSE)估计器：
<mrow> <msub> <mi>L</mi> <mrow> <mi>L</mi> <mi>M</mi> <mi>M</mi> <mi>S</mi> <mi>E</mi> </mrow> </msub> <mo>=</mo> <msub> <mi>L</mi> <mrow> <mi>T</mi> <mi>D</mi> <mi>C</mi> <mo>|</mo> <mi>&amp;mu;</mi> <mo>=</mo> <mn>1</mn> </mrow> </msub> <mo>=</mo> <msub> <mi>U</mi> <mrow> <mi>h</mi> <mn>1</mn> </mrow> </msub> <mfrac> <mrow> <msub> <mi>&amp;Lambda;</mi> <mrow> <mi>h</mi> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msubsup> <mi>&amp;sigma;</mi> <mrow> <mi>n</mi> <mi>o</mi> <mi>i</mi> <mi>s</mi> <mi>e</mi> </mrow> <mn>2</mn> </msubsup> <msub> <mi>E</mi> <mi>r</mi> </msub> </mrow> <msub> <mi>&amp;Lambda;</mi> <mrow> <mi>h</mi> <mn>1</mn> </mrow> </msub> </mfrac> <msubsup> <mi>U</mi> <mrow> <mi>h</mi> <mn>1</mn> </mrow> <mi>T</mi> </msubsup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>31</mn> <mo>)</mo> </mrow> </mrow>
当μ＝0，即为所谓最小二乘(LS)估计器：
<mrow> <msub> <mi>L</mi> <mrow> <mi>L</mi> <mi>S</mi> </mrow> </msub> <mo>=</mo> <msub> <mi>L</mi> <mrow> <mi>T</mi> <mi>D</mi> <mi>C</mi> <mo>|</mo> <mi>&amp;mu;</mi> <mo>=</mo> <mn>0</mn> </mrow> </msub> <mo>=</mo> <msub> <mi>U</mi> <mrow> <mi>h</mi> <mn>1</mn> </mrow> </msub> <msubsup> <mi>U</mi> <mrow> <mi>h</mi> <mn>1</mn> </mrow> <mi>T</mi> </msubsup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>32</mn> <mo>)</mo> </mrow> </mrow>
步骤6、求解纯净脉冲响应函数，根据式(30)可知纯净脉冲响应函数估计表示成
<mrow> <mover> <mi>h</mi> <mo>^</mo> </mover> <mo>=</mo> <msub> <mi>L</mi> <mrow> <mi>T</mi> <mi>D</mi> <mi>C</mi> </mrow> </msub> <mi>h</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>33</mn> <mo>)</mo> </mrow> </mrow>
步骤7、求解降噪后频响函数，在获得纯净脉冲响应函数估计后，再利用傅里叶变换即可得到降噪后频响函数；
步骤8、求解降噪结果，依据特征值分析结果，对比保留不同特征值个数对应的降噪效果，确定协方差特征值个数，估算其后验信噪比SNRdB，并取拉格朗日乘子μ，依据最优估计器LTDC求解得到降噪结果。
CN201710981385.8A 2017-10-20 2017-10-20 一种基于子空间特征值分解的实测频响函数降噪方法 Active CN107729845B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710981385.8A CN107729845B (zh) 2017-10-20 2017-10-20 一种基于子空间特征值分解的实测频响函数降噪方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710981385.8A CN107729845B (zh) 2017-10-20 2017-10-20 一种基于子空间特征值分解的实测频响函数降噪方法

Publications (2)

Publication Number Publication Date
CN107729845A true CN107729845A (zh) 2018-02-23
CN107729845B CN107729845B (zh) 2024-03-22

Family

ID=61212252

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710981385.8A Active CN107729845B (zh) 2017-10-20 2017-10-20 一种基于子空间特征值分解的实测频响函数降噪方法

Country Status (1)

Country Link
CN (1) CN107729845B (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108614804A (zh) * 2018-04-25 2018-10-02 中国人民解放军战略支援部队信息工程大学 基于信噪比检验的正则化卡尔曼滤波方法
CN110384490A (zh) * 2019-07-29 2019-10-29 杭州埃因霍温科技有限公司 基于相空间的bcg信号心率提取方法
CN111934716A (zh) * 2019-11-02 2020-11-13 广东石油化工学院 一种电力线通信信号滤波方法及系统
CN113204739A (zh) * 2021-05-24 2021-08-03 桂林电子科技大学 一种基于k-均值聚类的频响函数质量线优化方法
CN113686528A (zh) * 2021-07-28 2021-11-23 华南理工大学 一种结构-tld系统的子系统动力特性检测方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2006203890A (ja) * 2005-01-17 2006-08-03 Ntt Docomo Inc 周波数領域サブ空間チャネル推定装置及び方法、受信機及び信号受信方法
CN106646376A (zh) * 2016-12-05 2017-05-10 哈尔滨理工大学 基于加权修正参数的p范数噪声源定位识别方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2006203890A (ja) * 2005-01-17 2006-08-03 Ntt Docomo Inc 周波数領域サブ空間チャネル推定装置及び方法、受信機及び信号受信方法
CN106646376A (zh) * 2016-12-05 2017-05-10 哈尔滨理工大学 基于加权修正参数的p范数噪声源定位识别方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
孙鑫晖;张令弥;王彤;: "基于奇异值分解的频响函数降噪方法", 振动、测试与诊断, no. 03, 15 September 2009 (2009-09-15) *

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108614804A (zh) * 2018-04-25 2018-10-02 中国人民解放军战略支援部队信息工程大学 基于信噪比检验的正则化卡尔曼滤波方法
CN108614804B (zh) * 2018-04-25 2022-09-02 中国人民解放军战略支援部队信息工程大学 基于信噪比检验的正则化卡尔曼滤波方法
CN110384490A (zh) * 2019-07-29 2019-10-29 杭州埃因霍温科技有限公司 基于相空间的bcg信号心率提取方法
CN111934716A (zh) * 2019-11-02 2020-11-13 广东石油化工学院 一种电力线通信信号滤波方法及系统
CN111934716B (zh) * 2019-11-02 2021-09-17 广东石油化工学院 一种电力线通信信号滤波方法及系统
CN113204739A (zh) * 2021-05-24 2021-08-03 桂林电子科技大学 一种基于k-均值聚类的频响函数质量线优化方法
CN113686528A (zh) * 2021-07-28 2021-11-23 华南理工大学 一种结构-tld系统的子系统动力特性检测方法

Also Published As

Publication number Publication date
CN107729845B (zh) 2024-03-22

Similar Documents

Publication Publication Date Title
CN107729845A (zh) 一种基于子空间特征值分解的实测频响函数降噪方法
Cheng et al. A noise reduction method based on adaptive weighted symplectic geometry decomposition and its application in early gear fault diagnosis
Darong et al. A new incipient fault diagnosis method combining improved RLS and LMD algorithm for rolling bearings with strong background noise
CN102323518B (zh) 一种基于谱峭度的局部放电信号识别方法
CN105424366A (zh) 基于eemd自适应消噪的轴承故障诊断方法
CN110967599A (zh) 一种电能质量扰动检测与定位算法
CN107844627A (zh) 一种仅输出时变结构模态参数贝叶斯估计方法
Pan et al. A noise reduction method of symplectic singular mode decomposition based on Lagrange multiplier
CN103190898A (zh) 心磁信号噪声自适应滤波消除设计方法
CN110376575B (zh) 一种基于阻尼参数匹配随机共振的低频线谱检测方法
CN110749655B (zh) 一种针对比例阻尼结构的复模态辨识方法
CN113378661A (zh) 一种基于改进小波阈值和相关检测的直流电能信号去噪方法
CN106328155A (zh) 一种修正先验信噪比过估计的语音增强方法
CN106451498A (zh) 一种基于改进广义形态滤波的低频振荡模态辨识方法
CN103699513A (zh) 一种基于多尺度噪声调节的随机共振方法
CN107886078A (zh) 一种基于分层自适应阈值函数的小波阈值降噪方法
CN104427143B (zh) 残留回声检测方法及系统
Cui et al. Spectrum-based, full-band preprocessing, and two-dimensional separation of bearing and gear compound faults diagnosis
CN105738883A (zh) 一种部分均匀海杂波背景下的平滑广义似然比检测方法
Cheng et al. Enhanced symplectic characteristics mode decomposition method and its application in fault diagnosis of rolling bearing
CN110716532A (zh) 一种基于小波包能量与fft的水下机器人推进器弱故障辨识方法
CN113472390A (zh) 一种基于深度学习的跳频信号参数估计方法
CN101876585A (zh) 基于小波包估计噪声方差的ica收缩去噪方法
CN106980722B (zh) 一种脉冲响应中谐波成分的检测和去除方法
CN110865375B (zh) 一种水中目标检测方法

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