[发明专利]一种基于子空间特征值分解的实测频响函数降噪方法

摘要

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

主权项

一种基于子空间特征值分解的实测频响函数降噪方法，其特征在于：包括如下步骤：步骤1、频响函数的离散傅里叶反变换，子空间算法在时域内进行，根据傅里叶反变换将频响函数变为时域内的脉冲响应函数对于带噪的频响函数，h为脉冲响应函数，H与h均为离散观测序列，n为脉冲响应序列的长度,据离散傅里叶反变换hk=1nΣs=0n-1Hsej2πskn,0≤k≤n-1---(1)]]>将离散傅里叶反变换后的带噪脉冲响应函数写成由纯净的脉冲响应函数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)：Ns=Σk=0n-1nke-j2πks/n,nk=1nΣk=0n-1Nse+j2πks/n---(4)]]>将nk展开成实部与虚部nk=1nΣk=0n-1Nse+j2πks/n=1nΣk=0n-1Real(Ns)e+j2πks/n+1nΣk=0n-1Imag(Ns)e+j2πks/n---(5)]]>利用欧拉公式展开，经推导可证明nk也服从高斯分布，期望与方差分别为:E[nk]＝0  (6)Var[nk]=1n2nσ2+1n2nσ2=2nσ2=σnoise2---(7)]]>即对于频域带噪信号，若其噪声信号的实部与虚部独立，并服从均值为零的高斯分布，则经过傅里叶反变换后的脉冲响应函数同样服从高斯分布，期望与方差如上式；步骤2、构造线性滤波器L，将由纯净的脉冲响应函数h'k加上噪声nk构成的带噪脉冲响应函数成向量的形式：h＝h′+n  (8)上式中向量的长度均为n×1，假定h′的估计可通过h的线性滤波器L构成：h^=Lh---(9)]]>将残余误差写成由频响函数失真误差εh′和残余噪声误差εn两个部分构成其中，频响函数失真误差定义为：ϵh′=Δ(L-E)h′---(11)]]>残余噪声误差定义为：ϵn=ΔLn---(12)]]>由上述式(6)、(7)、以及aTa＝trace(aaT)推导出频响函数失真能量为：ϵh′2‾=E(ϵh′Tϵh′)=trace(E[ϵh′ϵh′T])=trace[(L-E)Rh′(L-E)T]---(13)]]>残余噪声能量为：ϵn2‾=E(ϵnTϵn)=trace(E[ϵnϵnT])=trace(LRnLT)=σnoise2trace(LLT)---(14)]]>其中，R*为协方差矩阵，trace为矩阵的迹，由上述推导可得到以下时域约束最优问题用于求解最优线性估计器：minLϵh′2‾]]>s.t.1nϵn2‾≤δ---(15)]]>其中,δ为正常数，即满足残余噪声能量小于δ时，频响函数失真最小，引入拉格朗日乘子μ，μ≥0，构造拉格朗日乘式：根据KKT条件:以及最终可得到最优估计器：LTDC=Rh′(Rh′+μRn)-1=Rh′(Rh′+μσnoise2Em)-1---(17)]]>步骤3、特征值分解，假定嵌入维数为m，延迟时间为τ，利用延迟坐标法重构脉冲响应函数的相空间矩阵，重构的相空间矩阵Y写成：对降噪分析而言，直接取延迟时间τ＝1，重构的相空间矩阵Y写成：其中，Y为Toeplitz矩阵，假设频响函数中存在加性噪声，则相空间矩阵Y可写成：Y＝X+N  (20)其中，Y、X和N分别表示原始带噪、纯净信号和噪声信号的脉冲响应函数的重构矩阵，带噪脉冲响应函数的协方差可写成：Rh＝Rh′+Rn  (21)采用下式对Rh进行估计：Rh≈1n-m+1YTY---(22)]]>步骤4、用Rh对Rh′进行估计，Rh为m×m矩阵,将纯净脉冲响应函数的协方差矩阵Rh′进行EVD分解：Rh′=Uh′Λh′Uh′T---(23)]]>将式(23)带入式(21)：Rh=Rh′+Rn=Uh′Λh′Uh′T+σnoise2Uh′EmUh′T=Uh′(Λh′+σnoise2Em)Uh′T---(24)]]>假定Rh′的秩为r，其中r＜m，则Λh′中的前r个对角元素不为零，即Rh′=[Uh′1Uh′2]Λh′1000Uh′1TUh′2T=Uh′1Λh′1Uh′1T---(25)]]>将式(25)带入(24)，则原始脉冲响应的协方差矩阵Rh可写成分块形式：Rh=[Uh′1Uh′2]Ah′1+σnoise2Er00σnoise2Em-rUh′1TUh′2T---(26)]]>直接对通过估计获得的Rh进行EVD分解，并按(26)进行分块：Rh=[Uh1Uh2]Λh100Λh2Uh1TUh2T=[Uh′1Uh′2]Λh100Λh2Uh′1TUh′2T---(27)]]>对比式(26)、(27)，可知：RhΛh1=Λh′1+σnoise2Er---(28)]]>将式(28)带入式(25)，可知：Rh′=Uh′1(Λh1-σnoise2Er)Uh′1T=Uh1(Λh1-σnoise2Er)Uh1T---(29)]]>步骤5、得出最优估计器，根据式(17)、(23)以及(29)，最终得最优估计器LTDCLTDC=Uh1Λh′1Λh′1+μσnoise2ErUh1T=Uh1Λh1-σnoise2ErΛh1+(μ-1)σnoise2ErUh1T---(30)]]>当μ＝1，式(30)即为所谓线性最小二乘误差(LMMSE)估计器：LLMMSE=LTDC|μ=1=Uh1Λh1-σnoise2ErΛh1Uh1T---(31)]]>当μ＝0，即为所谓最小二乘(LS)估计器：LLS=LTDC|μ=0=Uh1Uh1T---(32)]]>步骤6、求解纯净脉冲响应函数，根据式(30)可知纯净脉冲响应函数估计表示成h^=LTDCh---(33)]]>步骤7、求解降噪后频响函数，在获得纯净脉冲响应函数估计后，再利用傅里叶变换即可得到降噪后频响函数；步骤8、求解降噪结果，依据特征值分析结果，对比保留不同特征值个数对应的降噪效果，确定协方差特征值个数，估算其后验信噪比SNRdB，并取拉格朗日乘子μ，依据最优估计器LTDC求解得到降噪结果。