华中科技大学ADSP仿真作业(LMSLSLPHDMUSIC)

1、ADSP 实验报告实验报告院系：电子与通信学院作者：*学号：老师：日期：2015/11/20Huazhong University of Science and Technology 目目 录录ADSP 实验报告实验报告实验一实验一.11.1题目.11.2数据模型及参数选取.11.3算法模型.11.3.1自适应滤波器.11.3.2LMS算法简介.21.3.3LSL算法简介.21.4实验仿真.31.4.1仿真结果.31.4.2结果分析.51.5实验小结.5实验二实验二.62.1题目.62.2数据模型及参数选取.62.3算法模型.72.3.1PHD算法简介.72.3.2MUSIC算法简介.82.4

2、实验仿真.92.4.1仿真结果.92.4.2结果分析.102.5实验小结.11附录：源代码附录：源代码.12附录 A：实验一 .12附录 B：实验二.141实验一实验一1.1题目Implement LSL algorithm and LMS algorithm based on figure 3.30(P92) and figure 3.31(P93) .Model and parameters see page 91.课本 91 页模型及参数为：序列由 2 阶自回归模型产生：( )x n，是均值为 0、方差为 1 的高斯白)()2()1()(21nwnxanxanx( )w n噪声序列。用

3、LMS 算法和 LSL 算法来估计模型参数、121.558,0.81aa 1a，对两种算法的性能进行比较，并分析 LSL 算法中初始预测误差剩余对收2a敛性能的影响。1.2数据模型及参数选取1）高斯白噪声可以用中的函数产生均值为 0、方差为 1 的标matlabrandn准正态分布随机矩阵来实现，用模型产生信号，激励源是产生的高(2)ARrandn斯白噪声；2）信号点数 N=500，用 500 个信号来估计滤波器系数；3）依据题意，LMS 算法中选取，并取=0.1,1, 10来观察初始预0.005测误差剩余对 LSL 收敛性能的影响。1.3算法模型1.3.1自适应滤波器自适应滤波器自适应滤波器

4、由参数可调的数字滤波器和自适应算法两部分组成，如Error! Reference source not found.所示。输入信号通过参数可调的数字滤( )x n波器后产生输出信号，将其与参考信号进行比较，形成误差信号。( )y n( )d n( )e n通过某种自适应算法对滤波器参数进行调整，最终使的均方值最小。( )e n( )e n此时得到的是对的最佳估计。( )y n( )d n2自适应算法参数可调数字滤波器x(n)d(n)y(n)e(n)-+1.3-1 自适应滤波器原理1.3.2LMS 算法简介算法简介LMS 算法采用平方误差最小的原则代替最小均方误差最小的原则，使用LMS 算法来得

5、到权系数的公式如下：(1.3-1) nXnenwnw21其中：(1.3-2) nyndnenXnwnyT*式中，为时刻自适应滤波器的权矢量， w nn，为自适应滤波器的阶数；为时刻自 TNnwnwnwnw120,.,N X nn适应滤波器的参考输入矢量，由最近个信号采样值构成，N；是期望的输出值；为自适应 TNnxnxnxnX1,.1,( )d n( )e n滤波器的参考响应与输出相应之差，即误差信号；是控制自适应速度与稳定性的增益常数，又叫收敛因子或步长因子。1.3.3LSL 算法简介算法简介LSL 算法步骤如下：1）初始化0( )0,(0)0,(0)1,(0)(0)bffmmmme n这里

6、是前向和后向预测误差剩余的初始值。如果这个值没有给出，可以任意选择。2）迭代计算（按时间） ；1,2n 0020000( )( )( )( )( )(1)( )( )1bfbffe nenx ne nennxnn33）迭代计算（按阶）0,1,1mM11(1)( )( )(1)(1)bfmmmmmenennnn 11( )(1)( )( )(1)bffmmmmbmn enenenn11( )( )( )(1)( )fbbmmmmfmn enenenn211( )( )( )(1)ffmmmbmnnnn211( )( )(1)( )bbmmmfmnnnn21(1)(1)(1)(1)bmmmbmen

7、nnn同阶的嵌套着按时间进行迭代计算。1m各阶前向和后向反射系数可分别由以下二式计算：11( )( )( )fmmfmnknn11( )( )(1)bmmbmnknnM 是给定的滤波器的阶。采用 LSL 自适应算法，利用给出的数据可计算出格型滤波器的 1 阶和 2 阶前向和后向反射系数、和，然后由这些反射系数计算1( )fkn1( )bkn2( )fkn2( )bkn出估计、。经推导得1 a2 a 1112 ( )( )( )( )bfba nknkn kn （1.3-3） （1.3-4）22 ( )( )ba nkn1.4实验仿真1.4.1仿真结果仿真结果4图 1.4-1 AR(2)所产生信

8、号图 1.4-2 自适应权系数1 a、2 a的收敛轨迹 5图 1.4-3 初始预测误差剩余对 LSL 收敛性能的影响1.4.2结果分析结果分析1）由图 1.4-2 可看出，一是 LSL 算法和 LMS 算法计算得到的都收敛1 ( )a n于 1.558，都收敛于-0.81，但是 LSL 算法很明显地比 LMS 算法收敛得更2 ( )a n快，二是 LMS 算法得到的、的收敛轨迹有所波动，而 LSL 算法所得到收1 a2 a敛轨迹更加平稳；2）由图 1.4-3，在=0.1 和 1 两种情况下，很快就收敛到，而当1 ( )a n1a选为 10 时，收敛稍许慢了一些，但三种情况的收敛速度都要比 LM

9、S 算法快的多。1.5实验小结通过本次实验，对 LMS 算法和 LSL 算法的原理有了更深刻的理解，对影响 LMS 和 LSL 算法性能的一些参数也有了一定的认识，懂得了如何通过调节参数来获得想要的性能。6实验二题目给定，是方差为的复白噪声。信(2 *0.52)2 *0.54( )( )jnjnx neev n0,1,24n ( )v n2v噪比定义为211010log()vSNRdB1）用 PHD 方法估计两个频率，使用 M+1=3 阶自相关矩阵；2）MUSIC 方法估计频率，选取 N=12，M=2。2.2数据模型及参数选取1）依题意，可通过提供的函数可对复正弦信号( )x nmatlaba

10、wgn加高斯白噪声得到。(2 *0.52)2 *0.54( )jnjns nee2）PHD算法取样自相关点数_25Nsample 正弦信号个数2M 自相关矩阵阶数13NM 复白噪声功率 = 0.012v信噪比211010log20()vSNRdB待估计频率 0.5, 0.52f 3）MUSIC算法取样自相关点数_25Nsample 正弦信号个数2M 自相关矩阵阶数112NM 复白噪声功率 = 0.042v信噪比SNR =13.9794dB待估计频率 0.5, 0.52f 72.3算法模型2.3.1PHD 算法简介算法简介PHD（Pisarenko Harmonic Decomposition）

11、即谐波分解方法。在的特殊情况下，1NM是矩阵，是 M 阶对角阵，信号的自相关矩阵1 2MEeee(1)MM12mPdiag PPP可表示为SR HsREPE （2.3-1） 设是中的任一矢量，表示由映射x1MR1MR （2.3-2)sR xy得到的一组矢量，令y （2.3-3）HPE x于是 （2.3-4）1Mi iIye由此看出，的秩最多是 M，若可任意选择，则秩恰为 M。由于是矩阵，这SRi1M 就是要求他在中。利用式（2.3-3） ，只需证明是满秩的或=是线性独立的。MRHPE()HHPEEP但 （2.3-5）1Mii iiEPPe它不能是零，除非所有是零（假设所有是非零且所有的频率不同

12、） 。这就证明了的iPsR秩是且有唯一的零特征值。由于M （2.3-6）2xsvRRI与有相同的特征矢量，而的特征值是，所以的特征值是sRsR12 ,0M xR。现考虑对于最小特征值的特征矢量，有222212,vvMvvxR2v1Mv （2.3-7）211xMvMR vv由式（2.3-6） ，得到8 （2.3-8）1100sMHMR vEPE v左乘，得到1HMv （2.3-9）1111()()0HHHHHMMMMvEPE vE vP E v由于 P 是正定的，于是得出 （2.3-10）10HME v, （2.3-11）10HiMe v1,2,iM上式可写成展开形式, （2.3-12）10(1

13、)exp()0MMinvnjwn1,2,iM这就是说，如果给定，就能够求出多项式1Mv, （2.3-13）10(1),exp()MnMinvnzzjw1,2,iM位于单位圆上的零点，这些零点的角度就是正弦波的频率。2.3.2MUSIC 算法简介算法简介MUSIC（Multiple Signal Classification）即多信号分类算法，这种算法是以下列谱估计 （2.3-14）211( )|MUSICNHii MPfe v的峰作为正弦波频率的估计的。式中，表示频率， 表示信号矢量。理论上，当时，feiff，。由于存在估计误差，故在正弦波频率上或者附近将有一个峰，与式iee( )MUSICP

14、f （2.3-14）的 M 个最大峰对应的频率便是正弦波的频率估计。式（2.3-14）的最大值也可以利用信号子空间的特征矢量来进行计算。为此，注意到该式分母最小时最大，但分母可表示成 ( )MUSICPf （2.3-15）211 |NNHHHiiii Mi Me vevv e9该式利用了。式（2.3-14）的最大化等效于下式1 NHiiiIvv （2.3-16）21|NHii Me v的最大化。可见，无论是原先提出的噪声子空间特征矢量，还是式（2.3-15）的信号子空间矢量，都可以用来确定正弦波的频率估计。2.4实验仿真2.4.1仿真结果仿真结果PHD 算法中解方程结果（运行五次程序）：1f0

15、.45170.45910.45030.45200.45492f0.52550.53690.52470.52570.5298表 2.4-1 解方程求得的频率估计值图 2.4-1 作图法估计频率值101f0.43660.45330.44590.45050.46202f0.51500.52650.52600.52650.5403表 2.4-2 PHD 作图法捕捉到的频率值图 2.4-2 MUSIC 求估计频率值1f0.46940.47080.47350.47250.47122f0.51180.51230.51040.51500.5134表 2.4-3 MUSIC 方法捕捉到的频率值2.4.2结果分析

16、结果分析1）PHD 算法中，解方程的方法，求解方程（2.4-12）位于单位圆上的零点，并计算Z的角度，然后把变换到，再将角频率除以，就可以得到估计的频率值。Zww0,2 )w2f由表 2.1 可以看出，估计出来的两个频率值与实际值 0.5 和 0.52 存在一定的误差。PHD 的另一种估计频率的方法是，计算式（2.4-10）的值，在两个频率附近，值应该为0，但实际上由于估计误差，不会精确地等于 0，如果把计算的结果倒数，那么在估计的频率附近，将会出现尖峰，捕捉尖峰所对应的横轴值，就可以得到估计的频率值。由图 2.4-1、表2.4-4 与表 2.4-2 可以看出，两个峰值频点与实际都存在一定的误

17、差，但都在可以接受的范围11之内。2）图 2.4-2 为 MUSIC 方法估计频率，其方法与 PHD 的第二种方法类似，由表 2.4-5 可看出，该方法估计的频率值与标准值 0.5 和 0.52 之间也存在误差。2.5实验小结通过本次实验，对频率估计方法有了一定的了解，掌握了 PHD 和 MUSIC 算法估计频率方法。12附录：源代码附录：源代码附录 A：实验一clear;close all;% 由AR(2)模型生成输入信号% 初始化参数a1 = 1.558;a2 = -0.81; N = 500; %信号点数 %信号及白噪声信号序列的初始化x = zeros(1,N); %信号的初始化Noi

18、se = randn(1,N); %白噪声的初始化，均值为0，方差为1x(1) = Noise(1);x(2) = a1*x(1)+Noise(2);for n = 3:N x(n) = a1*x(n-1)+a2*x(n-2)+Noise(n);end% LMS算法估计a1、a2的收敛曲线 u = 0.005; %LMS算法下自适应增益常数初始化L = 2;LMS_w = zeros(L,n);for i = (L+1):n X = x(i-1:-1:(i-2); y(i) = X*LMS_w(:,i); %i时刻输出信号 e(i) = x(i)-y(i); LMS_w(:,(i+1) = L

19、MS_w(:,i)+2*u*e(i)*X; %i时刻滤波器的权值end LMS_a1 = LMS_w(1,1:n); %LMS算法中权系数a1的提取 LMS_a2 = LMS_w(2,1:n);% LSL算法估计a1、a2的收敛曲线 %初始化相关的参量M = 3;Delta(1:M,1) = 0; %前后向预测误差相关系数Sf(1:M,1) = 1; %初始前向预测误差剩余Sb(1:M,1) = 1; %初始后向预测误差剩余gama(1:M,1) = 1; %角参量%迭代计算for m = 1:M-113 for n = 2:N eb(1,n) = x(n); ef(1,n) = x(n);

20、Sb(1,n) = Sf(1,n-1)+x(n)*x(n); Sf(1,n) = Sf(1,n-1)+x(n)*x(n); gama(1,n) = 1; Delta(m+1,n) = Delta(m+1,n-1)+eb(m,n-1)*ef(m,n)/gama(m,n-1); ef(m+1,n) = ef(m,n)-Delta(m+1,n)*eb(m,n-1)/Sb(m,n-1); eb(m+1,n) = eb(m,n-1)-Delta(m+1,n)*ef(m,n)/Sf(m,n); Sf(m+1,n) = Sf(m,n)-Delta(m+1,n)*Delta(m+1,n)/Sb(m,n-1);

21、 Sb(m+1,n) = Sb(m,n-1)-Delta(m+1,n)*Delta(m+1,n)/Sb(m,n); gama(m+1,n-1) = gama(m,n-1)-eb(m,n-1)*eb(m,n-1)/Sb(m,n-1); kf(m+1,n) = Delta(m+1,n)/Sf(m,n); kb(m+1,n) = Delta(m+1,n)/Sb(m,n-1); endendkb(:,1) = 0;%初始kb，kf的值都是0，根据递推关系可以算出kf(:,1) = 0;for n=1:N LSL_a1(n)= kb(2,n)-kf(2,n)*kb(3,n); LSL_a2(n) = k

22、b(3,n);end% 绘图figure;plot(x);xlabel(n);ylabel(x(n);title(x(n) generated by AR model);figure();hold onplot(1:n,LMS_a1,g,1:n,LMS_a2,g);plot(1:n,LSL_a1,r,1:n,LSL_a2,r);plot(a1*ones(1,n),-);plot(a2*ones(1,n),-); text(0,1.65,a1=1.558);text(0,-0.85,a2=-0.81);text(100,1.2,LMS);text(400,1.7,LSL,delta=1);xla

23、bel(n);title(LSL与LMS算法估计信号模型参数的性能比较);axis(0 500 -1.5 2);hold off% 初始预测误差剩余对LSL收敛性能的影响D = 0.1,1,10; %初始预测误差剩余a1 = zeros(3,N); %不同初始值对应的a1的矩阵14%迭代计算for i = 1:3 Sf(1:M,1) = D(i); Sb(1:M,1) = D(i);for m = 1:M-1 for n = 2:N eb(1,n) = x(n); ef(1,n) = x(n); Sb(1,n) = Sf(1,n-1)+x(n)*x(n); Sf(1,n) = Sf(1,n-1

24、)+x(n)*x(n); gama(1,n) = 1; Delta(m+1,n) = Delta(m+1,n-1)+eb(m,n-1)*ef(m,n)/gama(m,n-1); ef(m+1,n) = ef(m,n)-Delta(m+1,n)*eb(m,n-1)/Sb(m,n-1); eb(m+1,n) = eb(m,n-1)-Delta(m+1,n)*ef(m,n)/Sf(m,n); Sf(m+1,n) = Sf(m,n)-Delta(m+1,n)*Delta(m+1,n)/Sb(m,n-1); Sb(m+1,n) = Sb(m,n-1)-Delta(m+1,n)*Delta(m+1,n)/

25、Sb(m,n); gama(m+1,n-1) = gama(m,n-1)-eb(m,n-1)*eb(m,n-1)/Sb(m,n-1); kf(m+1,n) = Delta(m+1,n)/Sf(m,n); kb(m+1,n) = Delta(m+1,n)/Sb(m,n-1); endendkb(:,1) = 0;%初始kb，kf的值都是0，根据递推关系可以算出kf(:,1) = 0;for n=1:N LSL_a1(n) = kb(2,n)-kf(2,n)*kb(3,n); a2(n) = kb(3,n);end a1(i,:) = LSL_a1;endfigure;hold onplot(a1

26、(1,:),g);plot(a1(2,:),b);plot(a1(3,:),r);axis(0 500 -0.5 2);xlabel(n);title(初始预测误差剩余delta对LSL收敛性能的影响);legend(delta=0.1,delta=1,delta=10);%附录 B：实验二1、PHD 算法clc; clear; % 初始值设定15N = 3; %自相关矩阵 Rx 的阶数 N_sample = 25; % 取样点数 M = 2; % 正弦信号的个数 sigma = 0.1; %白噪声的方差 SNR = 10 * log10(1/(sigma . 2); %求信噪比 % 产生输入

27、信号n = (0:N_sample-1); %数据取样点;s = exp(1i*2*pi*0.5*n)+exp(1i*(2*pi*0.52*n+0.25*pi); %正弦信号x = awgn(s,SNR); %向信号中加入白噪声% 实现 PHD 算法Rxx = xcorr(x); %求自相关函数 Rx = toeplitz(Rxx(N_sample:N_sample+N-1); %求自相关矩阵 V,D = eig(Rx); %求特征值和特征向量 V_min = V(:,1); %最小特征值对应的特征向量 %方法 1：解正交方程（2.3 12） ，直接求出频率值 Z = roots(V_min.); %方程（2.3 12）位于单位圆上的零点w = angle(Z); %零点的角度 for i = 1 : M if w(i) 0 w(i) = w(i) + 2 * pi % angle 函数求出的角度范围是-pi, pi), %当 w 时，需要加 2*pi,将范围调整到0, 2*pi)之间 endendf_PHD = w / (2*pi) %求得的频率值%方法 2：信号矢量的频率估计,通过方程（2.3 10）用作图法求估计频率值K = 1000; % 在 0-1 之间分 1000 步搜索 f = linspace(0, 1, K); e = zeros(N, K); for