ADSP仿真报告 华中科技大学.docx

1、ADSP仿真报告 华中科技大学专业：电磁场与微波技术学生成员：刘 航(M201571827) 完成时间：2020年1月14日图目录1 题目11.1 题目 Implement LSL algorithm and LMS algorithm based on figure 3.30(P92) and figure 3.31(P93). Model and parameters see page 91.P91模型如下：二阶自回归随机过程x(n)=a1x(n-1)+a2x(n-2)+w(n)，对应的滤波器参数a1=1.558，a2 =-0.81。分别用LMS算法和LSL算法预测信号x(n)，从而得到对

2、信号模型的两个参数值的估计值。1.2 算法模型1.2.1 自适应滤波原理自适应滤波器由参数可调的数字滤波器和自适应算法两部分组成，其原理如下图所示： 图 1.1自适应滤波器框图输入信号x(n)通过参数可调的数字滤波器后产生输出信号y(n)，将其与期望输出信号d(n)进行比较，形成误差信号e(n)。e(n)通过某种自适应算法对滤波器参数进行调整，最终使得e(n)的均方值最小。1.2.2 LMS算法简介在自适应算法中，最理想的情况是沿着性能曲面最陡的方向向下搜索曲面的最低点，每次迭代都需要知道性能曲面上某点的梯度值，而实际应用中梯度值只能根据观测数据进行估计。LMS（Least mean squa

3、re）最小均方算法，是一种更为简便的方法，它的核心思想是用平方误差代替均方误差，从而使得计算量大大减小，易于实现。LMS算法中，梯度用下式近似：可以得到：实际上，只是单个平方误差序列的梯度，而则是多个平方误差序列统计平均的梯度，所以LMS算法就是用前者作为后者的近似。得到LMS算法的基本关系式此时，全系数的调整过程是有“噪声”的，其估计路径不可能准确的沿着理想的最陡下降的路径。LMS算法中，信号基本关系总结如下：1.2.3 LSL算法简介基于最小均方误差(MMSE)准则的算法，如最陡下降法、LMS算法，其缺点是：收敛速度慢；对非平稳信号的适应性较差。为克服以上缺点，引入“最小二乘(LS)”准则

4、。最小二乘准则，是以误差的平方和最小作为最佳估计的一种误差准则。最小二乘滤波(LS滤波)可描述为：调整滤波器的权矢量，使得在每个时刻对所有已输入的信号而言，滤波器输出的误差平方和最小。它与LMS算法的区别是：LMS滤波是以集合平均为基础的，属于统计分析方法；LS滤波是一种瞬时分析方法, 即在每个时刻对所有已输入信号，都要重新评估其平方和，并通过调整权矢量使其达到最小。LSL是最小二乘格型算法（Least square lattice），用矢量空间法描述前向和后向线性预测误差，用预测误差滤波器的格型结构来递推更新误差参数。算法的具体推导比较复杂，详见教材的3.14节及3.15节。以下归纳LSL自

5、适应算法的计算过程与流程：(1) 初始化：；为一个小的正数，它是给定的前向和后向误差能量的初值, 若此值未预先给出,则可任意设定。(2) 迭代计算（按时间）；(3) 迭代计算（按阶数）；（这里注意同阶的 M+1嵌套着按时间的迭代运算）；。(4) 计算各阶计算前向和后向反射系数：；。1.3 本题模型对于本题，其模型框图如下：图 1.2 题目一模型根据照题中所给的2阶自回归模型，由单位方差高斯白噪声w(n)激励一个线性时不变全极点系统产生信号x(n)。预测误差为：y(n)=x(n)-d(n)，然后再按照某种准则控制预测误差，从而自适应的调节FIR滤波器的权系数，使之最终达到最优。本题用LMS算法以

6、及LSL算法对预测误差进行处理：LMS算法直接计出w1(n)及w2(n)；LSL算法计算出格型滤波器的1阶和2阶前向和后向反射系数，然后由这些系数计算出参数估计值，具体公式见教材3.244-3.277式。1.4 仿真过程及结果分析1.4.1 仿真过程(1) 初始化参数，由AR(2)模型生成x(n)。如下截图：其中信号点数取500，AR(2)模型由filter函数产生，filter函数具体用法参见matlab帮助。(2) LMS算法实现，程序如下截图：程序中for循环的公式见1.2.2节最后总结的信号基本关系的公式。矩阵w的第一行为参数a1的估计，第二行为参数a2的估计。(3)LSL算法实现，程

7、序截图如下：程序中要注意下标值，matlab数组下标从1开始，对于2阶情况，按公式迭代是从m=0,1，故程序中是从m=1,2迭代。下图中for循环实现的公式见1.2.3节总结的算法流程。计算得到的a1(n)和a2(n)为参数估计值。(3) 画图程序，比较两个算法的性能。程序截图如下：将两个算法得到的参数估计值的收敛轨迹画在同一张图中，用不同颜色表示出来。(4) 对LSL算法，前向和后向预测误差初始值选取不同值，画出收敛图。这部分直接在LSL算法程序中修改初始参数，选取=0.1,1.0,10三种情况，得到a1的估计值，然后在同一张图中画出来。1.4.2 结果分析(1) LMS和LSL算法性能比较

8、，程序图形如下：图 1.3 LSL和LMS算法估计信号模型参数的性能比较如图，4条颜色不同的线表示参数的收敛轨迹，黄色线为a1和a2的标准值。从图中可以看出，LSL算法和LMS算法计算得到的a1(n)都收敛于1.558，a2(n)都收敛于-0.81，但是LSL算法明显地比LMS算法收敛的更快。(2)如下图，对LSL算法，画出了对参数a1的估计a1(n)的收敛轨迹，分别取0.1,1.0,10三种情况。可以看出，在=0.1和1两种情况下（图中蓝色线和黑色线），a1(n)很快收敛于1.558，但是在轨迹的开始阶段，有一个下冲；=10时（红色线），收敛速度稍微慢一些，但开始阶段没有明显地下冲。图 1.

9、4 初始预测误差对LSL收敛性能的影响2 题目22.1 题目Implement PHD algorithm and MUSIC algorithm based on P164 4.25(4)&(5).4.25题目：计算机模拟实验，给定：x(n) = exp(j2 * 0.5n) + exp(j(2 * 0.52n + /4) + v(n), n = 0,1,2,24v(n)是方差为v2的复白噪声。信噪比定义为SNR = 10log10(1/v2)(dB)。(4)用PHD方法估计两个频率。使用M+1=3阶自相关矩阵，其元素是有偏自相关函数值。(5)用MUSIC方法估计频率，选取N=12，M=2。

10、计算式(4.155)时要采用更细的频率间隔。2.2 算法模型2.2.1 特征分解频率估计原理特征分解的主要思想是：把自相关矩阵的信息空间分成两个子空间，即信号子空间和噪声子空间；这两个子空间中的矢量的函数在正弦波频率上有尖锐的风，据此即可以估计正弦波的频率。总结如下：已知求正弦波频率估计。，其中Rs矩阵的秩为M，它有N-M个零特征值。所以Rx可以分解为两个子空间：噪声子空间（N-M维），它有N-M正交特征矢量，对应的特征值均为（即噪声功率）；信号子空间，有M个正交特征特征矢量，对应着M个大特征值。这里有一个重要性质：主特征矢量张成的子空间与由信号矢量张成的子空间是相同的。所以，信号矢量与噪声子

11、空间的所有矢量是正交的，即书上4.148式：这个性质是噪声子空间频率估计的基础。2.2.2 Music算法简介MUSIC算法就是基于上面讲的4.148式，定义下面这个谱函数（即书上4.155式）：式子中，f表示频率，e表示信号矢量，vi表示噪声子空间的特征矢量。即在处出现峰值，这些峰值所对应的M个f，就是所要估计的正弦信号频率。2.2.3 PHD算法简介在N=M+1的特殊情况下，是(M+1)xM矩阵，对Rx的最小特征值对应的特征矢量，有如下关系成立：所以，只要我们得到，就可以求出上式的解，这些解就是待估计的频率。2.3 仿真过程及结果分析2.3.1 仿真过程(1) 初始化参数，产生复白噪声和待

12、估计信号如图，根据题目要求，取样点数N_sample取25，信号个数M=2，由题目定义的信噪比产生复白噪声信号。(2) PHD算法实现A、 根据取样的观测信号x(n)，求出其N=M+1=3阶自相关矩阵。如下截图：这里使用xcorr函数和toeplitz函数。实际用的自相关矩阵的元素，用的是自相关函数的无偏估计（参见xcorr函数用法，第二个参数选unbiased）。用toeplitz函数，产生3阶自相关矩阵，保证得到的自相关矩阵具有理论自相关矩阵的性质（详见toeplitz函数的具体用法）。B、 求出Rx的特征值和特征向量，并找出最小特征值对应的特征向量C、 根据书上式4.152或者式4.15

13、4，求解方程，得到带估计频率值这里要注意，angle函数求出的角度范围是0, 2*pi),当求出的角频率w 0时，需要加上一个2*pi以调整到0, 2*pi)之间。其中roots函数具体用法参见matlab帮助。(3) MUSIC算法实现A、 第一步同PHD算法，求解自相关矩阵，并得到其特征矢量。其中，自相关矩阵的阶数N根据题目要求取N=12。B、 第二步，根据书上式4.155，用作图法，求出待估计频率代码中用到ginput函数，来认为捕捉图像的两个尖峰，然后输出了估计得到的两个频率。2.3.2 结果分析(1) 用PHD算法，运行得到频率的估计值：运行4次，得到的估计频率如下：，可以看出，PH

14、D算法能基本分辨出0.5和0.52两个频率，但误差较大。更改待估计信号的频率值，将两个信号的频率差距拉大，可以发现，得到的频率更准。如下（改成0.5和0.55两个频率）：,(2) 用MUSIC算法，捕捉图像的两个尖峰，得到频率的估计值：图 2.1 MUSIC算法得到的谱峰多运行几次，得到的频率值如下：，可以看出，MUSIC算法基本能够分辨出0.5和0.52两个频率。相比于PHD算法，其误差小一些，估计的更准一些。3 附录3.1 题目1源程序% 2阶自回归模型参数估计LSM&LSL算法by刘航 20151108clear;clf;clc;% 初始化参数%N是信号点数；N=500;na=1:N;%

15、二阶自回归信号模型a1=1.558;a2=-0.81;x=zeros(N,1);%Wa=randn(N,1);Wa=randn(1,N);%白噪声x=filter(1,1 -a1 -a2,Wa);% LMS算法过程L=2; %滤波器长度w=zeros(L,N); %LMS滤波器的系数mu=0.002; %mu是迭代步长；d=zeros(1,N);for i=L+1:N X=x(i-1) x(i-2); y(i)=X*w(:,i); e(i)=x(i)-y(i); w(:,i+1)=w(:,i)+2*mu*e(i)*X;end% LSL算法过程%初始化相关的参量M=3;Deta(1:M,1)=0

16、;gama(1:M,1)=1.0;sef(1:M,1)=1;seb(1:M,1)=1;%前向和后向预测误差初始值for m=1:M-1 %按阶数迭代（阶数为2,matlab下标从1开始，故M设为3） for n = 2 : N %包含按时间迭代 eb(1,n)=x(n); ef(1,n)=x(n); seb(1,n)=sef(1,n-1)+x(n)*x(n); sef(1,n)=sef(1,n-1)+x(n)*x(n); gama(1,n)=1; Deta(m+1,n)=Deta(m+1,n-1)+eb(m,n-1)*ef(m,n)/gama(m,n-1); ef(m+1,n)=ef(m,n)

17、-Deta(m+1,n)*eb(m,n-1)/seb(m,n-1); eb(m+1,n)=eb(m,n-1)-Deta(m+1,n)*ef(m,n)/sef(m,n); sef(m+1,n)=sef(m,n)-Deta(m+1,n)*Deta(m+1,n)/seb(m,n-1); seb(m+1,n)=seb(m,n-1)-Deta(m+1,n)*Deta(m+1,n)/seb(m,n); gama(m+1,n-1)=gama(m,n-1)-eb(m,n-1)*eb(m,n-1)/seb(m,n-1); kf(m+1,n)=Deta(m+1,n)/sef(m,n); kb(m+1,n)=Det

18、a(m+1,n)/seb(m,n-1); endendkb(:,1)=0;kf(:,1)=0;for n=1:N a1n(n)=kb(2,n)-kf(2,n)*kb(3,n); a2n(n)=kb(3,n);end% 画图程序figure(1),plot(na,w(1,na),k-,linewidth,2),hold on; plot(na,w(2,na),b-,linewidth,2),hold on;plot(na,a1n,g-,linewidth,2);hold on;plot(na,a2n,r-,linewidth,2);legend(表示LMS中a1,表示LMS中a2,表示LSL中a

19、1,表示LSL中a2,0);hold on;plot(1.558*ones(n),y-);plot(-0.81*ones(n),y-);text(0,1.558,1.558,color,r),text(0,-0.81,-0.81,color,r);title(LSL 与 LMS 性能参数比较) ;grid on;xlabel(n);ylabel(a1(n),a2(n);% LSL不同初始值收敛性能% %更改sef(1:M,1)=1;seb(1:M,1)=1的值为0.1和10(得到a1_1n,a1_2n)figure(2),plot(na,a1n,k-,linewidth,2),hold on;

20、plot(na,a1_1n,r-,linewidth,2),hold on;plot(na,a1_2n,b-,linewidth,2),hold on;legend(1.0,10,0.1,0);grid on;%3.2 题目2源程序% PHD算法by刘航20151108% 初始化参数，产生复白噪声噪声和待估计信号clc;clear all;N_sample = 25; % 取样自相关点数M = 2; % 正弦信号的个数sigma = 0.1; v_real = sigma * randn(N_sample, 1) / sqrt(2); % 白噪声的实部v_imag = sigma * rand

21、n(N_sample, 1) / sqrt(2); % 白噪声的虚部v = v_real + 1i*v_imag; % 复白噪声n = 1 : N_sample;n = n;s = exp(1i * 2 * pi * 0.5 * n) + exp(1i * (2 * pi * 0.55 * n + pi/4); x = s + v; % 观察信号% PHD算法实现%根据信号取样求出自相关序列，进而求出3阶相关矩阵N = 3; % 题(4)中自相关矩阵Rx的阶数x_corr,a = xcorr(x,unbiased); %求相关函数无偏估计Rx = (toeplitz(x_corr(N_samp

22、le : N_sample + N -1).; %求自相关矩阵% 求Rx的特征值和特征向量，并找出最小特征值对应的特征向量V,D = eig(Rx); % 求特征值和特征向量vmin = V(:,1); % 最小特征值对应的特征向量N=M+1% 解正交方程4.154，直接解出频率值Z = roots(vmin.); % 方程(4.152)求解，得到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*p

23、i)之间 endendf_PHD = w / (2*pi)% MUSCI算法by刘航20151109% 初始化clc;clear all;N = 12; % 题(5)中自相关矩阵Rx的阶数N_sample = 25; % 取样自相关点数M = 2; % 正弦信号的个数% 产生复白噪声噪声和信号sigma = 0.1; %白噪声功率SNR = 10 * log10(1/(sigma . 2); % 求信噪比v_real = sigma * randn(N_sample, 1) / sqrt(2); % 白噪声的实部v_imag = sigma * randn(N_sample, 1) / sqr

24、t(2); % 白噪声的虚部v = v_real + 1i*v_imag; % 复白噪声n = 1 : N_sample;n = n;s = exp(1i * 2 * pi * 0.5 * n) + exp(1i * (2 * pi * 0.52 * n + pi/4); x = s + v; % 观察信号% 数据自相关矩阵获得；x_corr = xcorr(x,unbiased); %相关序列的无偏估计unbiasedRx = (toeplitz(x_corr(N_sample : N_sample + N -1).; % 求自相关矩阵% 求Rx的特征值和特征向量V,D = eig(Rx);

25、 % 求特征值和特征向量vmin = V(:, 1:N-M); % 噪声子空间的N-M个特征向量% 通过方程(4.155)用作图法求估计频率值K = 10000; % 在0-1之间分10000步搜索f = linspace(0, 1, K);e = zeros(N, K);for n = 1 : N e(n,:) = exp(1i * (n-1) * 2 * pi * f);endy = e * vmin;y = (abs(y) . 2;sum_music = sum(y, 2);sum_music = 1 ./ sum_music;plot(f, sum_music);f_MUSIC, yy = ginput(2);% 捕捉两个尖峰，尖峰所对对应的横坐标即为带估计频率f_MUSIC