数字信号处理实验二

DSP Experiment II

1.Consider an AR process x(n) defined by the difference equation
$$x(n)=-a_{1}x(n-1)-a_{2}x(n-2)+v(n)$$
where v(n) is an additive white noise of zero mean and variance ${\sigma }_v^2$.The AR parameters $a_1$ and $a_2$ are both real valued:

$$\begin{gathered} a_1=0.1 \\ {a}_2=-0.8 \end{gathered}$$
a) Calculate the noise variance ${\sigma }_v^2$ such that the AR process x(n) has unit variance .Hence , generate different realization of the process x (n).

b) Given the input x (n), an LMS filter of length M = 2 is used to estimate the unknown AR parameters $a_1$ and $a_2$ . The step size $\delta$ is assigned the value 0.05. Compute and plot the ensemble average curve of $a_1$ and $a_2$ by averaging the value of parameters $a_1$ and $a_2$ over an ensemble of 100 different realization of the filter. Calculate the time constant according to the experiment results and compare with the corresponded theoretical value.

c) For one realization of the LMS filter, compute the prediction error
$$f(n)=x(n)-\tilde{x}(n)$$
And the two tap-weight errors
$${\epsilon }_1(n)=-a_1-h_1(n)$$
and
$${\epsilon }_2(n)=-a_2-h_2(n)$$
Using power spectral plots of $f(n),{\epsilon }_1(n),{\epsilon }_2(n)$ , show that $f(n)$ behaves as white noise, whereas ${\epsilon }_1(n),{\epsilon }_2(n)$ behave as low-pass process. Explain the reason for this phenomenon.

d) Compute the ensemble average learning curve of the LMS filter by averaging the square value of the prediction error f(n) over an ensemble of 100 different realization of the filter. Calculate the time constant and residual power according to the experiment results and compare with the corresponded theoretical values.

a)由于x(n)为一个AR过程，
$$R(N+1)\left[\begin{array}{c}1\\ -A_N\end{array}\right]=\left[\begin{array}{c}{\sigma }_v^2\\ 0\end{array}\right]$$

​ 其中
$$A_N=\left[\begin{array}{c}-a_1\\ -a_2\end{array}\right],R_{N+1}=\left[\begin{array}{ccc}r(0)& r(1)& r(2)\\ r(1)& r(0)& r(1)\\ r(2)& r(1)& r(0)\end{array}\right],r(0)=r_x^2-E^{2}x(n)=r_x^2=1$$

​ 将$a_1=0.1,a_2=-0.8$代入上式，可得$r(1)=-0.5,r(2)=0.85,{\sigma }_v^2=0.27$。将参数带入$x(n)$，实现图如下所示

b)由a)里面知道，
$$R_x=\left[\begin{array}{cc}r(0)& r(1)\\ r(1)& r(0)\end{array}\right]=\left[\begin{array}{cc}1& -0.5\\ -0.5& 1\end{array}\right]$$
​ 利用matlab可以得到 的两个特征值为 $\lambda_1=0.5 $ $\lambda_2=1.5 $。那么 $a_1$ 和 $a_2$ 对应的迭代曲线的理论的时间常数分别为${\tau }_1=1/\delta {\lambda }_1=13.3,{\tau }_2=1/\delta {\lambda }_2=40$ 。
a1=-0.1014 (第500点的值)
a2=1.0034 (第500点的值)

利用LMS算法
$$\{\begin{array}{l}e(n+1)=x(n+1)-W^T(n)X(n)\\ W(n+1)=W(n)+\delta e(n+1)X(n)\end{array}，H(n)=\left[\begin{array}{c}{h}_0(n)\\ h_1(n)\end{array}\right],X(n)=\left[\begin{array}{c}x(n)\\ x(n-1)\end{array}\right]$$

其中 $\delta =0.05$

100实验后取平均如图:

实际均值收敛时间常数可以根据时间常数的定义，从h0和h1曲线中取两个点计算时间常数。观测上图在 [0，10] 这一段曲线比较光滑，所以我们取第5点和第10点。$a(\infty )=[-0.07598,0.7414{]}^{\prime }$ 。因此根据公式
$$a(t1)-a(\infty )=e^{-\frac{t1-t2}{\acute{\tau }}}[a(t2)-a(\infty )]$$

得 $\acute{{\tau }_1}=18,\acute{{\tau }_2}=35.58$ 。

c)$f(n), \varepsilon _1(n), \varepsilon_2(n) $的功率谱见图。

从图中可以看出$a1(n),a2(n)$具有低通特性，$f(n)$具有白噪声特性。由于权系数最终会收敛到LMS的最优值，因此 $f(n)$ 会趋于 $v(n)$，具有白噪声特性。
而$E[{\alpha }_k(n)]=(1-\delta {\lambda }_k{)}^{n}E[{\alpha }_k(0)]$，系数将很快收敛到最优值a1和a2。由上式做DTFT
$$S(w)=\frac{E[{\alpha }_k(0)]}{1-(1-\delta {\lambda }_k)e^{-jw}},|S(w){|}^2=\frac{E^2[{\alpha }_k(0)]}{1+(1-\delta {\lambda }_k{)}^2-2(1-\delta {\lambda }_k)cos(w)}$$

上式具有低通特性， 即 $E[{\alpha }_k(n)]$ 具有低通特性，由$\alpha (n)=Q^T[H(n)-H_{opt}]$可知，$\alpha (n)$
是$h_i(n)$的线性组合，因此 $h_i(n)$ 具有低通特性，从而${\epsilon }_1(n)$和${\epsilon }_2(n)$具有低通特性

d)最陡下降法和LMS算法算法的均方误差曲线如图

由LMS算法性能分析可得，理论上均方收敛的时间为
$${\tau }_e=\frac{1}{2\delta \overline{\lambda }}=\frac{1}{2\times 0.05\times 1}=10$$

均方收敛后的误差为：$J(\infty ）=J_{min}(1+\frac{\delta }{2}N{\sigma }_x^2)$
将$J_{min}=R-H_{opt}^T{r}_{yx}$, $H_{opt}=R_{xx}^{-1}{r}_{yx}$, $R_{xx}=\left[\begin{array}{cc}{r}_x(0)& r_x(1)\\ r_x(1)& r_x(0)\end{array}\right]$, $r_{yx}=\left[\begin{array}{c}{r}_x(1)\\ r_x(2)\end{array}\right])$, N=2, ${\sigma }_x^2=r_x(0)$
代入上式，得到理论上最小均方差$J_{min}=0.2700$均方收敛误差为：$J(\infty )=0.2842$

​ 实验中，$J(\infty )$的波动还是很大，取收敛的几百个点取平均，近似求得：$J(\infty)=0.28$91在上图的曲线中取，$J(15)=0.5197, J(47)=0.2955, J(\infty)=0.2833$, 由公式 $J(15)-J(\infty )=e^{-\frac{15-47}{\tau }}[J(47)-J(\infty )]$，可得$\tau \approx 10$

​ 综上所述，实验所得的均方收敛时间和均方收敛误差和理论值基本吻合。由于超量误差$J_{ex}\approx J(\infty )-J_{min}$ ，而 $J_{min}=0.2700$，因此可得理论上的超量误差为 $J_{ex}\approx 0.0142$ ，实验中的超量误差为 $J_{ex}=0.0191$

---

$x_1(n)$ 自相关函数为：
Rx1(k)=E[x1(n)x1(n−k)]=E{sin(0.05πn+φ)sin[0.05π(n−k)+φ]}=−0.5×E{cos[0.05π(2n−k)+2φ]−cos(0.05πk)}=0.5cos(0.05πk) \begin{aligned} R_{x_1}(k)&=E[x_1(n)x_1(n-k)]\\ &=E\{sin(0.05\pi n+\varphi)sin[0.05\pi(n-k)+\varphi] \}\\ &=-0.5\times E\{cos[0.05\pi(2n-k)+2\varphi] -cos(0.05\pi k) \}\\ &=0.5cos(0.05\pi k) \end{aligned} Rx1​​(k)​=E[x1​(n)x1​(n−k)]=E{sin(0.05πn+φ)sin[0.05π(n−k)+φ]}=−0.5×E{cos[0.05π(2n−k)+2φ]−cos(0.05πk)}=0.5cos(0.05πk)​
​ $x_2(n)$ 自相关函数为：
$$R_{x_2}(m)=\{\begin{array}{ll}6& m=0\\ 4& m=1\\ 1& m=2\\ 0& m\ge 3\end{array}$$

​ 如图所示：

由上图可见，当延时 $D>3$ 时，就可以保证LMS滤波器参考信号中的宽带信号成分与
输入信号中的宽带信号成分不相关，而两者中的窄带信号成分仍然保持一定的相关性，因
此通过LMS滤波后能够有效地起到谱线增强的作用。

​ 考虑到噪声功率比信号功率强太多，取滤波器阶数高，步长短。 滤波器长度选N=400 , $\sigma$ 取0.000002。

如上图所示：可以较好的恢复原有的信号。