维纳滤波

介绍这两种算法之前先来简单介绍下维纳滤波的问题

x(n)x\left( n \right)x(n)和y(n)y\left( n \right)y(n)是零均值的平稳离散信号，并且已知它们的二阶矩，采用最小均方误差准则（MMSE, Minimum Mean-Squared Error）可以设计最优线性滤波器。
均方误差JJJ为
J=E[e2(n)]=E{[y(n)−y~(n)]2}J=E\left[ {{e}^{2}}\left( n \right) \right]=E\left\{ {{\left[ y\left( n \right)-\tilde{y}\left( n \right) \right]}^{2}} \right\}J=E[e2(n)]=E{[y(n)−y~(n)]2}
假设滤波器的冲激响应为h(n)h\left( n \right)h(n)，为了方便，将其记作hn{{h}_{n}}hn，则滤波器输出的结果为
y~(n)=hn∗x(n)=∑ihix(n−i)\tilde{y}\left( n \right)={{h}_{n}}*x\left( n \right)=\sum\limits_{i}^{{}}{{{h}_{i}}x\left( n-i \right)}y~(n)=hn∗x(n)=i∑hix(n−i)
而
e(n)=y(n)−y~(n)=y(n)−∑ihix(n−i)e\left( n \right)=y\left( n \right)-\tilde{y}\left( n \right)=y\left( n \right)-\sum\limits_{i}^{{}}{{{h}_{i}}x\left( n-i \right)}e(n)=y(n)−y~(n)=y(n)−i∑hix(n−i)
那么可以对均方误差关于冲激响应求导
∂J∂hj=2E[e(n)∂e(n)∂hj]=−2E[e(n)x(n−j)]=0\frac{\partial J}{\partial {{h}_{j}}}=2E\left[ e\left( n \right)\frac{\partial e\left( n \right)}{\partial {{h}_{j}}} \right]=-2E\left[ e\left( n \right)x\left( n-j \right) \right]=0∂hj∂J=2E[e(n)∂hj∂e(n)]=−2E[e(n)x(n−j)]=0
因此可以得到
E[e(n)x(n−j)]=0E\left[ e\left( n \right)x\left( n-j \right) \right]=0E[e(n)x(n−j)]=0
同样也可以得到E[e(n)y~(n)]=0E\left[ e\left( n \right)\tilde{y}\left( n \right) \right]=0E[e(n)y~(n)]=0
接下来，可以定义
rc(j)=E[y(n)x(n−j)]{{r}_{c}}\left( j \right)=E\left[ y\left( n \right)x\left( n-j \right) \right]rc(j)=E[y(n)x(n−j)]
r(j)=E[x(n)x(n−j)]r\left( j \right)=E\left[ x\left( n \right)x\left( n-j \right) \right]r(j)=E[x(n)x(n−j)]
可以得到
rc(j)=∑ihir(j−i){{r}_{c}}\left( j \right)=\sum\limits_{i}^{{}}{{{h}_{i}}r\left( j-i \right)}rc(j)=i∑hir(j−i)
这就是著名的Weiner-Hopf方程，线性最优滤波（维纳滤波）等价于将参考信号分为两路正交分量（误差信号和滤波器输出信号），并且误差信号和输入信号正交，滤波器输出信号与输入信号不正交。
考虑NNN阶FIR维纳滤波器的解，系统的输出可以表示为
X(n)=[x(n),x(n−1),⋯,x(n−N+1)]T\mathbf{X}\left( n \right)={{\left[ x\left( n \right),x\left( n-1 \right),\cdots ,x\left( n-N+1 \right) \right]}^{T}}X(n)=[x(n),x(n−1),⋯,x(n−N+1)]T
滤波器的系数
H=[h0,h1,⋯,hN−1]T\mathbf{H}={{\left[ {{h}_{0}},{{h}_{1}},\cdots ,{{h}_{N-1}} \right]}^{T}}H=[h0,h1,⋯,hN−1]T
则滤波器的输出可以表示为
y~(n)=HTX(n)=XT(n)H\tilde{y}\left( n \right)={{\mathbf{H}}^{T}}\mathbf{X}\left( n \right)={{\mathbf{X}}^{T}}\left( n \right)\mathbf{H}y~(n)=HTX(n)=XT(n)H
由正交性条件E[e(n)x(n−j)]=0E\left[ e\left( n \right)x\left( n-j \right) \right]=0E[e(n)x(n−j)]=0可以得到
E[e(n)X(n)]=0E\left[ e\left( n \right)\mathbf{X}\left( n \right) \right]=0E[e(n)X(n)]=0
E{[y(n)−HTX(n)]X(n)}=0E\left\{ \left[ y\left( n \right)-{{\mathbf{H}}^{T}}\mathbf{X}\left( n \right) \right]\mathbf{X}\left( n \right) \right\}\text{=0}E{[y(n)−HTX(n)]X(n)}=0
即
E[y(n)X(n)]−E[X(n)XT(n)]H=0E\left[ y\left( n \right)\mathbf{X}\left( n \right) \right]-E\left[ \mathbf{X}\left( n \right){{\mathbf{X}}^{T}}\left( n \right) \right]\mathbf{H}=0E[y(n)X(n)]−E[X(n)XT(n)]H=0
那么最优解可以表示为
Hopt=Rxx−1ryx{{\mathbf{H}}_{\text{opt}}}=\mathbf{R}_{xx}^{-1}{{\mathbf{r}}_{yx}}Hopt=Rxx−1ryx
其中
Rxx=[r(0)r(1)⋯r(N−1)r(1)r(0)⋯r(N−2)⋮⋮⋱⋮r(N−1)r(N−2)⋯r(0)]{{\mathbf{R}}_{xx}}=\left[ \begin{matrix} r\left( 0 \right) & r\left( 1 \right) & \cdots & r\left( N-1 \right) \\ r\left( 1 \right) & r\left( 0 \right) & \cdots & r\left( N-2 \right) \\ \vdots & \vdots & \ddots & \vdots \\ r\left( N-1 \right) & r\left( N-2 \right) & \cdots & r\left( 0 \right) \\ \end{matrix} \right]Rxx=⎣⎢⎢⎢⎡r(0)r(1)⋮r(N−1)r(1)r(0)⋮r(N−2)⋯⋯⋱⋯r(N−1)r(N−2)⋮r(0)⎦⎥⎥⎥⎤
ryx=[rc(0),rc(1),⋯,rc(N−1)]{{\mathbf{r}}_{yx}}=\left[ {{r}_{c}}\left( 0 \right),{{r}_{c}}\left( 1 \right),\cdots ,{{r}_{c}}\left( N-1 \right) \right]ryx=[rc(0),rc(1),⋯,rc(N−1)]
rc(p)=E[y(n)x(n−p)]{{r}_{c}}\left( p \right)=E\left[ y\left( n \right)x\left( n-p \right) \right]rc(p)=E[y(n)x(n−p)]
均方误差JJJ可以表示为
J=E[e2(n)]=E{[y(n)−y~(n)]2}=E[y2(n)]−2E[y(n)HTX(n)]+E[HTX(n)XT(n)H]=E[y2(n)]−2ryxTH+HTRxxH\begin{aligned} & J=E\left[ {{e}^{2}}\left( n \right) \right]=E\left\{ {{\left[ y\left( n \right)-\tilde{y}\left( n \right) \right]}^{2}} \right\} \\ & =E\left[ {{y}^{2}}\left( n \right) \right]-2E\left[ y\left( n \right){{\mathbf{H}}^{T}}\mathbf{X}\left( n \right) \right]+E\left[ {{\mathbf{H}}^{T}}\mathbf{X}\left( n \right){{\mathbf{X}}^{T}}\left( n \right)\mathbf{H} \right] \\ & =E\left[ {{y}^{2}}\left( n \right) \right]-2\mathbf{r}_{yx}^{T}\mathbf{H}+{{\mathbf{H}}^{T}}{{\mathbf{R}}_{xx}}\mathbf{H} \end{aligned} J=E[e2(n)]=E{[y(n)−y~(n)]2}=E[y2(n)]−2E[y(n)HTX(n)]+E[HTX(n)XT(n)H]=E[y2(n)]−2ryxTH+HTRxxH
最小化均方误差为
Jmin⁡=E[y2(n)]−2ryxTHopt+HoptTRxxHopt=E[y2(n)]−HoptTRxxHopt=E[y2(n)]−HoptTryx\begin{aligned} & {{J}_{\min }}=E\left[ {{y}^{2}}\left( n \right) \right]-2\mathbf{r}_{yx}^{T}{{\mathbf{H}}_{\text{opt}}}+\mathbf{H}_{\text{opt}}^{T}{{\mathbf{R}}_{xx}}{{\mathbf{H}}_{\text{opt}}} \\ & =E\left[ {{y}^{2}}\left( n \right) \right]-\mathbf{H}_{\text{opt}}^{T}{{\mathbf{R}}_{xx}}{{\mathbf{H}}_{\text{opt}}} \\ & =E\left[ {{y}^{2}}\left( n \right) \right]-\mathbf{H}_{\text{opt}}^{T}{{\mathbf{r}}_{yx}} \\ \end{aligned} Jmin=E[y2(n)]−2ryxTHopt+HoptTRxxHopt=E[y2(n)]−HoptTRxxHopt=E[y2(n)]−HoptTryx
误差性能曲面可以表示为
J=E[y2(n)]−2ryxTH+HTRxxH=Jmin⁡+(H−Hopt)TRxx(H−Hopt)\begin{aligned} & J=E\left[ {{y}^{2}}\left( n \right) \right]-2\mathbf{r}_{yx}^{T}\mathbf{H}+{{\mathbf{H}}^{T}}{{\mathbf{R}}_{xx}}\mathbf{H} \\ & ={{J}_{\min }}+{{\left( \mathbf{H}-{{\mathbf{H}}_{\text{opt}}} \right)}^{T}}{{\mathbf{R}}_{xx}}\left( \mathbf{H}-{{\mathbf{H}}_{\text{opt}}} \right) \end{aligned}J=E[y2(n)]−2ryxTH+HTRxxH=Jmin+(H−Hopt)TRxx(H−Hopt)

最陡梯度下降算法

对于维纳滤波问题来说，就是为了求最优的滤波器系数，同时有
J(Hopt)≤J(H)J\left( {{\mathbf{H}}_{\text{opt}}} \right)\le J\left( \mathbf{H} \right)J(Hopt)≤J(H)
那么，可以通过构造迭代算法来更新滤波器的系数H(0),H(1),⋯,H(k+1)\mathbf{H}\left( 0 \right),\mathbf{H}\left( 1 \right),\cdots ,\mathbf{H}\left( k+1 \right)H(0),H(1),⋯,H(k+1)，使得均方误差随着滤波器系数矢量的更新而下降，即
J[H(k)]>J[H(k+1)]J\left[ \mathbf{H}\left( k \right) \right]>J\left[ \mathbf{H}\left( k+1 \right) \right]J[H(k)]>J[H(k+1)]
可以采用沿着均方误差的最陡下降方向，即沿着均方误差的梯度方向的相反方向来更新滤波器的系数。
均方误差相对于系数矢量的梯度
VG(n)=∂E[e2(n)]∂H∣H=H(n)=∂J(H)∂H∣H=H(n)=−2E[e(n)X(n)]=2RxxH(n)−2ryx\begin{aligned} & {{\mathbf{V}}_{G}}\left( n \right)=\frac{\partial E\left[ {{e}^{2}}\left( n \right) \right]}{\partial \mathbf{H}}\left| _{\mathbf{H}=\mathbf{H}\left( n \right)} \right.=\frac{\partial J\left( \mathbf{H} \right)}{\partial \mathbf{H}}\left| _{\mathbf{H}=\mathbf{H}\left( n \right)} \right. \\ & =-2E\left[ e\left( n \right)\mathbf{X}\left( n \right) \right]=2{{\mathbf{R}}_{xx}}\mathbf{H}\left( n \right)-2{{\mathbf{r}}_{yx}} \end{aligned} VG(n)=∂H∂E[e2(n)]∣∣H=H(n)=∂H∂J(H)∣∣H=H(n)=−2E[e(n)X(n)]=2RxxH(n)−2ryx
那么滤波器系数的更新可以表示为
H(n+1)=H(n)−12δVG(n)= H(n)+δ[ryx−RxxH(n)]\mathbf{H}\left( n+1 \right)=\mathbf{H}\left( n \right)-\frac{1}{2}\delta {{\mathbf{V}}_{G}}\left( n \right)\text{=}\ \mathbf{H}\left( n \right)\text{+}\delta \left[ {{\mathbf{r}}_{yx}}-{{\mathbf{R}}_{xx}}\mathbf{H}\left( n \right) \right]H(n+1)=H(n)−21δVG(n)= H(n)+δ[ryx−RxxH(n)]
其中，δ>0\delta >0δ>0是进行迭代的步长。
从上面的过程可以看出，最陡梯度下降方法取决于滤波器的初值H(0)\mathbf{H}\left( 0 \right)H(0)，梯度VG{{\mathbf{V}}_{G}}VG，迭代步长δ\deltaδ，在迭代步长足够小的时候其收敛于维纳最优解。

最小均方误差(Least Mean Square,LMS)算法

虽说维纳滤波是最优解，但是其需要一定的先验知识，那么在缺少先验知识的时候该方法就不再适用。LMS自适应是一种重要的梯度自适应算法，它是用瞬时值代替估计值的一种方法，LMS具体的过程为
步骤一：给定初值H(0)\mathbf{H}\left( 0 \right)H(0)
步骤二：迭代
{e(n+1)=y(n+1)−HT(n)X(n+1)H(n+1)=H(n)+δe(n+1)X(n+1)\left\{ \begin{aligned} & e\left( n+1 \right)=y\left( n+1 \right)-{{\mathbf{H}}^{\text{T}}}\left( n \right)\mathbf{X}\left( n+1 \right) \\ & \mathbf{H}\left( n+1 \right)\text{=}\mathbf{H}\left( n \right)+\delta e\left( n+1 \right)\mathbf{X}\left( n+1 \right) \\ \end{aligned} \right. {e(n+1)=y(n+1)−HT(n)X(n+1)H(n+1)=H(n)+δe(n+1)X(n+1)

二阶FIR的LMS自适应滤波器

下图是用一个二阶FIR的LMS自适应滤波器消除正弦干扰的一个方案。

数值分析

最陡下降算法

H(n+1)=H(n)−12δVG(n)n=0,1,2,⋯\mathbf{H}\left( n+1 \right)=\mathbf{H}\left( n \right)-\frac{1}{2}\delta {{\mathbf{V}}_{G}}\left( n \right)\ \ \ \ \ \ \ n=0,1,2,\cdotsH(n+1)=H(n)−21δVG(n) n=0,1,2,⋯
式中，δ=0.4\delta =0.4δ=0.4表示迭代的步长，H(0)=[3,−4]T\mathbf{H}\left( 0 \right)={{[3,-4]}^{\text{T}}}H(0)=[3,−4]T
梯度计算
VG(n)=∂E[e2(n)]∂H∣H=H(n)=∂J(H)∂H∣H=H(n)=−2E[e(n)X(n)]=2RxxH(n)−2ryx\begin{aligned} & {{\mathbf{V}}_{G}}\left( n \right)=\frac{\partial E\left[ {{e}^{2}}\left( n \right) \right]}{\partial \mathbf{H}}\left| _{\mathbf{H}=\mathbf{H}\left( n \right)} \right.=\frac{\partial J\left( \mathbf{H} \right)}{\partial \mathbf{H}}\left| _{\mathbf{H}=\mathbf{H}\left( n \right)} \right. \\ & =-2E\left[ e\left( n \right)\mathbf{X}\left( n \right) \right]=2{{\mathbf{R}}_{xx}}\mathbf{H}\left( n \right)-2{{\mathbf{r}}_{yx}} \end{aligned} VG(n)=∂H∂E[e2(n)]∣∣H=H(n)=∂H∂J(H)∣∣H=H(n)=−2E[e(n)X(n)]=2RxxH(n)−2ryx
式中
Rxx(k)=116∑i=015(2sin⁡2πi16)(2sin⁡2π(i−k)16)=cos⁡2πk16{{\mathbf{R}}_{xx}}(k)=\frac{1}{16}\sum\limits_{i=0}^{15}{\left( \sqrt{2}\sin \frac{2\pi i}{16} \right)}\left( \sqrt{2}\sin \frac{2\pi (i-k)}{16} \right)=\cos \frac{2\pi k}{16}Rxx(k)=161i=0∑15(2sin162πi)(2sin162π(i−k))=cos162πk
ryx(k)=116∑i=015[sin⁡(2πi16+π10)](2sin⁡2π(i−k)16)=12cos⁡(2πk16+π10){{\mathbf{r}}_{yx}}(k)=\frac{1}{16}\sum\limits_{i=0}^{15}{\left[ \sin \left( \frac{2\pi i}{16}+\frac{\pi }{10} \right) \right]}\left( \sqrt{2}\sin \frac{2\pi (i-k)}{16} \right)=\frac{1}{\sqrt{2}}\cos \left( \frac{2\pi k}{16}+\frac{\pi }{10} \right)ryx(k)=161i=0∑15[sin(162πi+10π)](2sin162π(i−k))=21cos(162πk+10π)
则可以求出
VG(0)=2[Rxx(0)Rxx(1)Rxx(1)Rxx(0)][H(0)H(1)]−2[ryx(0)ryx(1)]=2[10.92390.92391][3−4]−2[0.67250.5377]=[−2.7326−3.5320]\begin{aligned} & {{\mathbf{V}}_{G}}\left( 0 \right)=2\left[ \begin{matrix} {{\mathbf{R}}_{xx}}\left( 0 \right) & {{\mathbf{R}}_{xx}}\left( 1 \right) \\ {{\mathbf{R}}_{xx}}\left( 1 \right) & {{\mathbf{R}}_{xx}}\left( 0 \right) \\ \end{matrix} \right]\left[ \begin{aligned} & \mathbf{H}\left( 0 \right) \\ & \mathbf{H}\left( 1 \right) \\ \end{aligned} \right]-2\left[ \begin{aligned} & {{r}_{yx}}\left( 0 \right) \\ & {{r}_{yx}}\left( 1 \right) \\ \end{aligned} \right] \\ & =2\left[ \begin{matrix} 1 & 0.9239 \\ 0.9239 & 1 \\ \end{matrix} \right]\left[ \begin{matrix} 3 \\ -4 \\ \end{matrix} \right]-2\left[ \begin{aligned} & 0.6725 \\ & 0.5377 \\ \end{aligned} \right] \\ & =\left[ \begin{aligned} & -2.7326 \\ & -3.5320 \\ \end{aligned} \right] \end{aligned}VG(0)=2[Rxx(0)Rxx(1)Rxx(1)Rxx(0)][H(0)H(1)]−2[ryx(0)ryx(1)]=2[10.92390.92391][3−4]−2[0.67250.5377]=[−2.7326−3.5320]
同时可以计算出
E[y2(n)]=ryy(0)=rss(0)+rN0N0(0)=0.05+0.5=0.55E\left[ {{y}^{2}}\left( n \right) \right]={{\mathbf{r}}_{yy}}\left( 0 \right)={{\mathbf{r}}_{ss}}\left( 0 \right)+{{\mathbf{r}}_{{{N}_{0}}{{N}_{0}}}}\left( 0 \right)=0.05+0.5=0.55E[y2(n)]=ryy(0)=rss(0)+rN0N0(0)=0.05+0.5=0.55
得到最优滤波器的系数计算公式
Hopt=Rxx−1ryx=[h0∗,h1∗]T{{\mathbf{H}}_{\text{opt}}}=\mathbf{R}_{xx}^{-1}{{\mathbf{r}}_{yx}}\text{=}{{\left[ {{h}_{0}}^{*},h_{1}^{*} \right]}^{\text{T}}}Hopt=Rxx−1ryx=[h0∗,h1∗]T
将各个初值带入计算可以得到
[h0∗h1∗]=[Rxx(0)Rxx(1)Rxx(1)Rxx(0)]−1[ryx(0)ryx(1)]=[1.200−0.571]\left[ \begin{aligned} & h_{0}^{*} \\ & h_{1}^{*} \\ \end{aligned} \right]={{\left[ \begin{matrix} {{\mathbf{R}}_{xx}}\left( 0 \right) & {{\mathbf{R}}_{xx}}\left( 1 \right) \\ {{\mathbf{R}}_{xx}}\left( 1 \right) & {{\mathbf{R}}_{xx}}\left( 0 \right) \\ \end{matrix} \right]}^{-1}}\left[ \begin{aligned} & {{\mathbf{r}}_{yx}}\left( 0 \right) \\ & {{\mathbf{r}}_{yx}}\left( 1 \right) \\ \end{aligned} \right]=\left[ \begin{matrix} 1.200 \\ -0.571 \\ \end{matrix} \right] [h0∗h1∗]=[Rxx(0)Rxx(1)Rxx(1)Rxx(0)]−1[ryx(0)ryx(1)]=[1.200−0.571]
误差函数为
J(n)=E[y2(n)]−2ryxTH+HTRxxH=0.55+h02+h12+2h0h1cos⁡π8−2h0cos⁡π10−2h1cos⁡9π40\begin{aligned} & J(n)=E\left[ {{y}^{2}}(n) \right]-2\mathbf{r}_{yx}^{T}\mathbf{H}+{{\mathbf{H}}^{T}}{{\mathbf{R}}_{xx}}\mathbf{H} \\ & =0.55+h_{0}^{2}+h_{1}^{2}+2{{h}_{0}}{{h}_{1}}\cos \frac{\pi }{8}-\sqrt{2}{{h}_{0}}\cos \frac{\pi }{10}-\sqrt{2}{{h}_{1}}\cos \frac{9\pi }{40} \end{aligned}J(n)=E[y2(n)]−2ryxTH+HTRxxH=0.55+h02+h12+2h0h1cos8π−2h0cos10π−2h1cos409π
此时可以求出最小误差
Jmin⁡=E[y2(n)]−HoptTRxxHopt=ryy(0)−HoptTryx=0.05{{J}_{\min }}=E\left[ {{y}^{2}}\left( n \right) \right]-\mathbf{H}_{\text{opt}}^{T}{{\mathbf{R}}_{xx}}{{\mathbf{H}}_{opt}}={{\mathbf{r}}_{yy}}\left( 0 \right)-\mathbf{H}_{\text{opt}}^{T}{{\mathbf{r}}_{yx}}=0.05Jmin=E[y2(n)]−HoptTRxxHopt=ryy(0)−HoptTryx=0.05
可以根据上面的分析画出误差性能曲面如下：

从图中可以看出，误差性能曲面是存在最小值的，这与上面的分析也一致。为了更加方便的观察，也可以给出其二维平面图，即等高线

很明显的看出中间就是最小值的地方。
为了更好地说明这两种算法的区别，下面给出两种算法的搜索过程

从图中可以看出，最陡梯度下降方法能够很好地接近最终的结果，但是LMS算法却在最终结果附近来回波动，同时在搜索的时候LMS算法的波动范围也比较大，可以通过增加实验次数来减小这种现象，100次实验结果如下

从图中可以看出，经过多次实验后LMS算法有了明显的改善，其搜索过程逐渐向最终的解靠近，这也与上面分析一致。
代码如下：

clear;
close all;
clc;
N=1000;                                                             %输入信号的点数
n=0:1:N-1;
iter=100;                                                             %实验次数
H_g1=zeros(N,iter);
H_g2=zeros(N,iter);
H0_a=zeros(N-1,iter);                                           %存放h0
H1_a=zeros(N-1,iter);                                           %存放h1
for j=1:itersn=sqrt(0.05)*randn(1,N);                                    %产生均值为0，方差为0.05的高斯白噪声xn=sqrt(2)*sin(2*pi*n/16);                                    %xn的表达形式yn=sn+sin(2*pi*n/16+pi/10);                               %yn的表达形式delta=0.4;                                                            %定义迭代步长%% 最陡梯度下降delta=0.4;                                                             %定义迭代步长H_g=zeros(2,N);                                                    %初始化滤波器矩阵大小H_g(:,1)=[3 -4];                                                      %滤波器初始值V_G=[-2.7326;-3.5320];                                          %初始化V_GRxx=[1 0.9239;0.9239 1];                                        %初始化RxxRyx=[0.6725;0.5377];                                              %初始化Ryxfor i=1:NV_G=2.*Rxx*H_g(:,i)-2.*Ryx;                                 %计算V_G的迭代值H_g(:,i+1)=H_g(:,i)-1/2.*delta.*V_G;                     %计算滤波器的更新值H_g1(i,j)=H_g(1,i);H_g2(i,j)=H_g(2,i);end%% LMS算法H_L=zeros(2,N);                                                    %初始化滤波器矩阵大小H_L(:,1)=[3;-4];                                                     %初始化h0,h1的值e=zeros(1,N);                                                      %存放误差的矩阵for i=1:N-1e(1,i)=yn(i+1)-H_L(:,i).'*[xn(i+1);xn(i)];               %计算每次的误差H_L(:,i+1)=H_L(:,i)+delta.*e(1,i).*[xn(i+1);xn(i)];H0_a(i,j)=H_L(1,i);H1_a(i,j)=H_L(2,i);end
end
%% 平均值
for i=1:N-1H0_g(i)=sum(H_g1(i,:))./iter;H1_g(i)=sum(H_g2(i,:))./iter;H0(i)=sum(H0_a(i,:))./iter;H1(i)=sum(H1_a(i,:))/iter;
end
[h0,h1]=meshgrid(-2:0.1:4,-4:0.1:2);                      %形成坐标格点矩阵
J_n=0.55+h0.^2+h1.^2+2*h0.*h1.*cos(pi/8)-sqrt(2).*h0.*cos(pi/10)-sqrt(2).*h1.*cos(9*pi/40);            %误差函数表达式
figure(1);
v=0:0.2:2;                                                            %画出高度为v的等值线
contour(h0,h1,J_n,v,'linewidth',1.5);
set(gca,'xtick',-2:1:3);                                            %设置坐标
set(gca,'ytick',-4:1:2);
xlabel('h0','fontsize',20);ylabel('h1','position',[-2.5 -1],'fontsize',20,'rotation',360);
hold on;
plot(H0_g,H1_g,'k','linewidth',3)
plot(H0,H1,'r','linewidth',3);
xlim([-2 4]);
h0_i=1.200;h1_i=-0.571;
line([-2,4],[h1_i,h1_i],'linewidth',2);                         %画出点（h0,h1）
line([h0_i,h0_i],[-4,2],'linewidth',2);
plot(h0_i,h1_i,'b*')
l=legend('误差等高线','最陡梯度下降算法','LMS算法');
set(l,'fontsize',10);

最陡梯度下降算法和LMS算法原理介绍及MATLAB实现相关推荐

浅谈迪杰斯特拉（Dijkstra）算法和A*算法原理及实现
写在前面最近我在学习一门名叫<智能自主机器人及系统>的课程,虽然跟过去所学的<机器人学>在部分内容上有所重复,但该课程的应用性更强.对于不同的机器人,如差速轮式车.四轮车.四 ...
Dijkstra算法和Floyd算法详解（MATLAB代码）
一.Dijkstra算法 1.算法简介 Dijkstra算法是由E.W.Dijkstra于1959年提出,又叫迪杰斯特拉算法,它应用了贪心算法模式,是目前公认的最好的求解最短路径的方法.算法解决的是有 ...
图像处理：TDLMS算法原理介绍及MATLAB实现
一.TDLMS介绍 1.1 算法原理二维最小均方(two-dimensional least mean square, TDLMS)滤波算法由最小均方误差(least mean square err ...
最差性能最优波束形成算法原理介绍及MATLAB实现
原理介绍标准的Capon波束形成是最小化阵列输出功率,并且对期望信号无失真接收来达到最大化输出信干噪比的目的.但是当真实导向矢量未知且存在一个不确定集中时,为了对真实导向矢量设置无失真响应约束,最差 ...
【数据挖掘实验】关联规则——CARMA算法和AprioriAll算法
一.实验项目名称: 关联规则--CARMA算法和AprioriAll算法二.实验目的与要求: 在软件方面:会用Clementine软件进行序列关联规则分析. 在理论方面:CARMA算法和Aprior ...
使用Apriori算法和FP-growth算法进行关联分析
目录 1. 关联分析 2. Apriori原理 3. 使用Apriori算法来发现频繁集 4. 使用FP-growth算法来高效发现频繁项集 5. 示例:从新闻网站点击流中挖掘新闻报道扩展阅读系列 ...
数据结构与算法之美笔记——基础篇（下）：图、字符串匹配算法（BF 算法和 RK 算法、BM 算法和 KMP 算法、Trie 树和 AC 自动机）
图如何存储微博.微信等社交网络中的好友关系?图.实际上,涉及图的算法有很多,也非常复杂,比如图的搜索.最短路径.最小生成树.二分图等等.我们今天聚焦在图存储这一方面,后面会分好几节来依次讲解图相关的 ...
关联规则挖掘算法: Aprior算法和Fpgrowth算法
关联规则挖掘的目的是挖掘不同物品(item)之前的相关性,啤酒和尿布的故事就是一个典型的成功例子.关联规则挖掘思想简单高效,在广告推荐领域也有较多的应用,主要用于推荐模型落地前的流量探索以及构建规 ...
BF算法和KMP算法
给定两个字符串S和T,在主串S中查找子串T的过程称为串匹配(string matching,也称模式匹配),T称为模式.这里将介绍处理串匹配问题的两种算法,BF算法和KMP算法. BF算法 (暴力匹配 ...

最陡梯度下降算法和LMS算法原理介绍及MATLAB实现