现代信号处理1

现代信号处理Assignment4姓名：饶康班级：研1305学号：2013020141题目1：考虑如下图权值线性组合器，输入端引入随机信号kr，其平均功率为20.01kEr；假设信号随机抽样相互独立，取16N。编程实现：（1）画出LMS算法性能曲面等值线，要求等值线权值间隔不超过1，标明坐标值、均方误差值和性能最小点位置及最小均方误差值，分别对应初始权值010,0.100ww和014,0.0510ww绘出加权值收敛轨迹，迭代次数不小于100次；（2）计算0.05和0.10时学习曲线的时间常数，绘出学习曲线并在学习曲线中观测时间常数，与理论计算值比较；（3）计算0.05和0.10时的失调并比较；（4）分析比较μ的大小对自适应滤波的影响。要求：写出实验报告：包括原理、方法和结果，并附源代码（加必要的注释）和仿真数据结果。解：（1）LMS自适应滤波器是使滤波器的输出信号与期望响应之间的误差的均方值为最小，因此称为最小均方（LMS）自适应滤波器。如题图，LMS自适应滤波器是有线性组合器构成，输入为x(k)和x(k-1)，输出y(k)是这些输入加权后的线性组合，即y(k)=𝑤0𝑥(𝑘)+𝑤1𝑥(𝑘−1)(1)定义权向量为W=[w0w1]T，输入信号向量为X=[xkxk-1]，题图中dk代表期望响应，定义误差kε+-w1w0z-1krNπk2sinNπkdk2cos2XkYk信号为𝜀𝑘=𝑑(𝑘)−𝑦(𝑘)(2)写成向量形式为𝜀𝑘=𝑑(𝑘)−𝑋𝑊(3)误差平方为𝜀𝑘2=𝑑2(𝑘)−2𝑑(𝑘)𝑋𝑊+𝑊𝑇𝑋𝑇𝑋𝑊(4)上式两边去数学期望后，得均方误差为E{𝜀𝑘2}=𝐸{𝑑2(𝑘)}−2𝐸{𝑑(𝑘)𝑋}𝑊+𝑊𝑇𝐸{𝑋𝑇𝑋}𝑊(5)互相关函数行向量定义为𝑅𝑑𝑥𝑇:𝑅𝑑𝑥𝑇=E{𝑑(𝑘)𝑋}(6)和自相关函数矩阵𝑅𝑥𝑥=E{𝑋𝑇𝑋}(7)则均方误差可表述为E{𝜀𝑘2}=𝐸{𝑑2(𝑘)}−2𝑅𝑑𝑥𝑇𝑊+𝑊𝑇𝑅𝑥𝑥𝑊(8)这表明，均方误差是权系数向量W的二次函数，它是一个中间向上凹的抛物型曲面，是具有唯一最小值的函数。调节权系数使均方误差为最小，相当于沿抛物形曲面下降找最小值。将(8)式对权系数W求导，得到均方误差函数的梯度为∇(k)=−2𝑅𝑑𝑥+2𝑅𝑥𝑥𝑊(9)另∇(k)=0，即可求出最佳权系数向量𝑊𝑜𝑝𝑡=𝑅𝑥𝑥−1𝑅𝑑𝑥(10)得到最优权值即可求最小均方误差，带入（8）式可得𝐸{𝜀𝑘2}𝑚𝑖𝑛=𝐸{𝑑2(𝑘)}−𝑅𝑑𝑥𝑇𝑊𝑜𝑝𝑡(11)以上计算最优权值的过程需要知道自相关和互相关，还要对矩阵求逆，还是比较麻烦的，而WidrowandHoffLMS算法给出了一种用迭代方法求最优权值近似值的方法，权值迭代公式如下W(k+1)=W(k)+2μ𝜀𝑘𝑋𝑇(12)其中μ是一个控制收敛速度与稳定性的常数，称之为收敛因子（或增益常数）。有了迭代公式就可以在给定初始权值以及收敛因子的条件下求出最优权值的近似值。以上算法使用matlab实现，初始权值选择W=[00]，μ=0.1，进行迭代运算200次，得到权值迭代轨迹，然后改变初始权值W=[4−10]，μ=0.05，再次进行迭代并再次得到权值轨迹。接下来计算题中系统的性能曲面，输入为单频正弦波，求其自相关可等效成整周期内在时间上求平均，即E{𝑥𝑘𝑥𝑘}=1𝑁∑(𝑠𝑖𝑛2𝜋𝑘𝑁)2+𝐸{𝑟𝑘2}𝑁−1𝑘=0=0.51E{𝑥𝑘𝑥𝑘−1}=1𝑁∑𝑠𝑖𝑛2𝜋𝑘𝑁𝑠𝑖𝑛2𝜋(𝑘−1)𝑁+𝐸{𝑟𝑘𝑟𝑘−1}=0.5cos(𝜋/8)𝑁−1𝑘=0再计算出期望信号和输入信号的互相关值，E{𝑑𝑘𝑥𝑘}=1𝑁∑2𝑐𝑜𝑠2𝜋𝑘𝑁𝑠𝑖𝑛2𝜋𝑘𝑁=0𝑁−1𝑘=0E{𝑑𝑘𝑥𝑘−1}=1𝑁∑2𝑐𝑜𝑠2𝜋𝑘𝑁𝑠𝑖𝑛2𝜋(𝑘−1)𝑁=−sin(𝜋/8)𝑁−1𝑘=0则可写出相关矩阵如下𝑅𝑥𝑥=[0.510.5cos(𝜋/8)0.5cos(𝜋/8)0.51]𝑅𝑑𝑥=[0−sin(𝜋/8)]最优权值为𝑊∗=𝑅𝑥𝑥−1𝑅𝑑𝑥[𝑤0∗𝑤1∗]=[0.510.5cos(𝜋/8)0.5cos(𝜋/8)0.51][0−sin(𝜋/8)]=[3.784−4.178]接着计算最小均方误差为𝜉𝑚𝑖𝑛=𝐸{𝑑2(𝑘)}−𝑅𝑑𝑥𝑇𝑊∗=0.4011那么性能函数为MSE=𝛏=𝜉𝑚𝑖𝑛+(𝑊−𝑊∗)𝑇𝑅𝑥𝑥(𝑊−𝑊∗)用matlab计算该函数并画出性能曲面的等值线，在图上标注好均方误差、性能最小点以及最小均方误差，然后将上述迭代的权值再画入图中得到如图1结果图1性能曲面等值线图（2）学习曲线的时间常数定义为𝜏𝑚𝑠𝑒=14𝜇(13)当收敛因子分别为0.05和0.1时计算时间常数，见图3。而学习曲线定义为均方误差MSE对应迭代次数的曲线，采用初始权值为W=[00]，收敛因子μ分别为0.1和0.05时计算学习曲线MSE(k)，并分别绘于图2(a)和图2(b)中如下图2(a)μ=0.05时学习曲线图2(b)μ=0.1时学习曲线从图上可以看出μ越小学习曲线收敛的越慢，也就是说学习曲线的时间常数越小，这与理论计算结果是一致的。（3）失调定义为M=μtr(𝑅𝑥𝑥)(14)其中tr(𝑅𝑥𝑥)是自相关矩阵的迹，即矩阵主对角线元素的和。分别对应收敛因子μ为0.05和0.1时计算出失调如图3图3时间常数和失调的计算（4）μ的取值会影响自适应滤波的收敛性，μ的取值不能太大。当满足0μ1𝜆𝑚𝑎𝑥(15)时才会收敛，其中𝜆𝑚𝑎𝑥为自相关矩阵的最大自相关值。而μ的大小还会影响收敛速度，较大的μ值会使学习曲线的时间常数更大，从而收敛更快，如图2中(a)图收敛速度比(b)图快。但是μ过大会形成震动性的过渡特性，就像图2中(a)图的收敛μ值较小，学习曲线比较平滑，(b)图的μ值较大，就会出现振动的特性。题目2：利用matlab设计一个二阶自适应噪声对消器，对含高斯白噪声信道干扰的正弦信号进行滤波，信噪比为10dB。画出有噪信号、降噪后信号和噪声收敛结果。解：首先设计噪声对消系统模型如图4图4自适应噪声对消系统模型噪声V0(n)、V1(n)、信源s(n)和y(n)是不相关的，但是噪声V0(n)和V1(n)之间是相关的。自适应滤波系统的误差为ε=s(n)+𝑉0(𝑛)−𝑦(n)(16)取平方的期望后E[𝜀2]=E[𝑠(𝑛)2]+E[(𝑉0(𝑛)−𝑦(𝑛))2](17)当调节滤波器使均方误差最小时，信号功率E[s2]将不受影响，即𝐸𝑚𝑖𝑛[𝜀2]=𝐸[𝑠(𝑛)2]+𝐸𝑚𝑖𝑛[(𝑉0(𝑛)−𝑦(𝑛))2](18)所以当均方误差最小时E[(𝑉0(𝑛)−𝑦(𝑛))2]也最小，此时便达到减弱或对消噪声的效果，此时的误差ε就是降噪后的信号。自适应滤波部分使用LMS算法，权值迭代的公式如题1中的迭代公式，经过matlab编程得到自适应噪声对消的效果如图5图5自适应噪声对消y(n)信源-d(n)噪声V0(n)∑∑∑AF∑+++噪声V1(n)s(n)e(n)附录一：题目1的matlab源代码%LMS算法closeallclearallclc%设置参数，N为采样个数，u为增益常数N=16;u=[0.10.05];%设置迭代次数kK=400;k=1:K;%计算出的自行关阵Rx=[0.510.5*cos(pi/8);0.5*cos(pi/8)0.51];Rdx=[0;-sin(pi/8)];rk=0+0.1*randn(1,K);%平均功率为0.01的随机信号（选用高斯噪声）dk=2*cos(2*pi*k/N);%期望信号xk=sin(2*pi*k/N)+rk;%输入的带噪信号%设置起始权值wk{1}(:,1)=[0;0];wk{2}(:,1)=[4;-10];%求理论最佳权值wk_l=(Rx^-1)*Rdx;%性能曲面等值线的等值e=[2015105210.80.60.50.42];symsw0symsw1v=[w0;w1]-wk_l;fori=1:2%用LMS算法迭代求最佳权值forj=2:Kxkk=[xk(j),xk(j-1)];%输入信号yk=xkk*wk{i}(:,j-1);%输出信号err=dk(j)-yk;%误差wk{i}(:,j)=wk{i}(:,j-1)+2*u(i)*err*xkk';%权值迭代endend%性能曲面最小值（最小均方误差）Emin=mean(dk.^2)-Rdx'*wk_l;figure(1)form=1:10z=Emin+v'*Rx*v-e(m);set(ezplot(z,[-28-100]),'color','blue')holdon;end%画迭代时权值的变化plot(wk{1}(1,:),wk{1}(2,:),'r');holdonplot(wk{2}(1,:),wk{2}(2,:),'r');holdon%标注最佳权值的位置plot(wk_l(1),wk_l(2),'*');title('性能曲面等值线')%计算学习曲线fori=1:2wk{1}(:,1)=[0;0];wk{2}(:,1)=[0;0];forj=2:Kxk=[sin(2*pi*j/N)sin(2*pi*(j-1)/N)]+rk(j);%输入信号yk=xk*wk{i}(:,j-1);%输出信号err=dk(j)-yk;%误差wk{i}(:,j)=wk{i}(:,j-1)+2*u(i)*err*xk';%权值迭代endfigureforj=1:Kmse(j)=Emin+(wk{i}(:,j)-wk_l)'*Rx*(wk{i}(:,j)-wk_l);endplot(mse)title('LMS算法学习曲线')xlabel('迭代次数k')ylabel('学习曲线MSE')ifi==1legend('u=0.1')elselegend('u=0.05')endend%计算时间常数t=1./(4.*u);%计算失调M=u.*trace(Rx);disp(['增益常数为0.1时，时间常数为',num2str(t(1)),',失调为',num2str(M(1))])disp(['增益常数为0.05时，时间常数为',num2str(t(2)),',失调为',num2str(M(2))])附录二：题目2的matlab源代码%用LMS算法设计自适应滤波器closeallclearallclcT=10000;t=0:0.001:0.001*(T-1);xs=sin(2*pi*t);xn=randn(1,T);xn=xn-mean(xn);power=1/length(xs)*sum(xs.*xs);var=power/10;xn=sqrt(var)/std(xn)*xn;%噪声，作为系统输入x=xs+xn;%带噪声信号e1=sum(xs.^2);e2=sum((x-xs).^2);snr=10*log10(e1/e2);d=xs+xn;%对于自适应对消器，用xs+xn作为期望信号%设计自适应滤波器N=2;%滤波器阶数w=zeros(N,1);%初始化滤波器权值rx_max=max(eig(xn*xn'));%输入信号相关矩阵的最大特征值u=abs(randn(1)*(1/rx_max));%收敛因子y=zeros(1,T);fork=2:Ty(k)=[xn(k)xn(k-1)]*w;e(k)=d(k)-y(k);w=w+2*u*e(k)*[xn(k)xn(k-1