ADSP仿真报告 华中科技大学.docx
《ADSP仿真报告 华中科技大学.docx》由会员分享，可在线阅读，更多相关《ADSP仿真报告 华中科技大学.docx（21页珍藏版）》请在冰点文库上搜索。
ADSP仿真报告华中科技大学
专业：
电磁场与微波技术
学生成员：
刘航（M201571827）
完成时间：
2020年1月14日
图目录
1
题目1
1.1题目
ImplementLSLalgorithmandLMSalgorithmbasedonfigure3.30（P92）andfigure3.31（P93）.Modelandparametersseepage91.
P91模型如下：
二阶自回归随机过程x（n）=a1x（n-1）+a2x（n-2）+w（n），对应的滤波器参数a1=1.558，a2=-0.81。
分别用LMS算法和LSL算法预测信号x（n），从而得到对信号模型的两个参数值的估计值。
1.2算法模型
1.2.1自适应滤波原理
自适应滤波器由参数可调的数字滤波器和自适应算法两部分组成，其原理如下图所示：
图1.1自适应滤波器框图
输入信号x（n）通过参数可调的数字滤波器后产生输出信号y（n），将其与期望输出信号d（n）进行比较，形成误差信号e（n）。
e（n）通过某种自适应算法对滤波器参数进行调整，最终使得e（n）的均方值最小。
1.2.2LMS算法简介
在自适应算法中，最理想的情况是沿着性能曲面最陡的方向向下搜索曲面的最低点，每次迭代都需要知道性能曲面上某点的梯度值，而实际应用中梯度值只能根据观测数据进行估计。
LMS（Leastmeansquare）最小均方算法，是一种更为简便的方法，它的核心思想是用平方误差代替均方误差，从而使得计算量大大减小，易于实现。
LMS算法中，梯度用下式近似：
可以得到：
实际上，
只是单个平方误差序列的梯度，而
则是多个平方误差序列统计平均的梯度，所以LMS算法就是用前者作为后者的近似。
得到LMS算法的基本关系式
此时，全系数的调整过程是有“噪声”的，其估计路径不可能准确的沿着理想的最陡下降的路径。
LMS算法中，信号基本关系总结如下：
1.2.3LSL算法简介
基于最小均方误差（MMSE）准则的算法，如最陡下降法、LMS算法，其缺点是：
收敛速度慢；对非平稳信号的适应性较差。
为克服以上缺点，引入“最小二乘（LS）”准则。
最小二乘准则，是以误差的平方和最小作为最佳估计的一种误差准则。
最小二乘滤波（LS滤波）可描述为：
调整滤波器的权矢量，使得在每个时刻对所有已输入的信号而言，滤波器输出的误差平方和最小。
它与LMS算法的区别是：
LMS滤波是以集合平均为基础的，属于统计分析方法；LS滤波是一种瞬时分析方法,即在每个时刻对所有已输入信号，都要重新评估其平方和，并通过调整权矢量使其达到最小。
LSL是最小二乘格型算法（Leastsquarelattice），用矢量空间法描述前向和后向线性预测误差，用预测误差滤波器的格型结构来递推更新误差参数。
算法的具体推导比较复杂，详见教材的3.14节及3.15节。
以下归纳LSL自适应算法的计算过程与流程：
（1）初始化：
；
；
；
；
为一个小的正数，它是给定的前向和后向误差能量的初值,若此值未预先给出,则可任意设定。
（2）迭代计算（按时间
）
；
；
；
（3）迭代计算（按阶数
）
；（这里注意同阶的M+1嵌套着按时间的迭代运算）
；
；
；
；
。
（4）计算各阶计算前向和后向反射系数：
；
。
1.3本题模型
对于本题，其模型框图如下：
图1.2题目一模型
根据照题中所给的2阶自回归模型，由单位方差高斯白噪声w（n）激励一个线性时不变全极点系统产生信号x（n）。
预测误差为：
y（n）=x（n）-d（n），然后再按照某种准则控制预测误差，从而自适应的调节FIR滤波器的权系数，使之最终达到最优。
本题用LMS算法以及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算法实现，程序截图如下：
程序中要注意下标值，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算法性能比较，程序图形如下：
图1.3LSL和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.4初始预测误差
对LSL收敛性能的影响
2题目2
2.1题目
ImplementPHDalgorithmandMUSICalgorithmbasedonP1644.25（4）&（5）.
4.25题目：
计算机模拟实验，给定：
x（n）=exp（j2π*0.5n）+exp（j（2π*0.52n+π/4））+v（n）,n=0,1,2,…,24
v（n）是方差为δv2的复白噪声。
信噪比定义为SNR=10log10（1/δv2）（dB）。
（4）用PHD方法估计两个频率。
使用M+1=3阶自相关矩阵，其元素是有偏自相关函数值。
（5）用MUSIC方法估计频率，选取N=12，M=2。
计算式（4.155）时要采用更细的频率间隔。
2.2算法模型
2.2.1特征分解频率估计原理
特征分解的主要思想是：
把自相关矩阵的信息空间分成两个子空间，即信号子空间和噪声子空间；这两个子空间中的矢量的函数在正弦波频率上有尖锐的风，据此即可以估计正弦波的频率。
总结如下：
已知
求正弦波频率估计。
，其中Rs矩阵的秩为M，它有N-M个零特征值。
所以Rx可以分解为两个子空间：
噪声子空间（N-M维），它有N-M正交特征矢量，对应的特征值均为
（即噪声功率）；信号子空间，有M个正交特征特征矢量，对应着M个大特征值
。
这里有一个重要性质：
主特征矢量
张成的子空间与由信号矢量
张成的子空间是相同的。
所以，信号矢量
与噪声子空间的所有矢量是正交的，即书上4.148式：
这个性质是噪声子空间频率估计的基础。
2.2.2Music算法简介
MUSIC算法就是基于上面讲的4.148式，定义下面这个谱函数（即书上4.155式）：
式子中，f表示频率，e表示信号矢量，vi表示噪声子空间的特征矢量。
即
在
处出现峰值，这些峰值所对应的M个f，就是所要估计的正弦信号频率。
2.2.3PHD算法简介
在N=M+1的特殊情况下，
是（M+1）xM矩阵，对Rx的最小特征值
对应的特征矢量
，有如下关系成立：
所以，只要我们得到
，就可以求出上式的解，这些解就是待估计的频率。
2.3仿真过程及结果分析
2.3.1仿真过程
（1）初始化参数，产生复白噪声和待估计信号
如图，根据题目要求，取样点数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.154，求解方程，得到带估计频率值
这里要注意，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次，得到的估计频率如下：
，
，
，
可以看出，PHD算法能基本分辨出0.5和0.52两个频率，但误差较大。
更改待估计信号的频率值，将两个信号的频率差距拉大，可以发现，得到的频率更准。
如下（改成0.5和0.55两个频率）：
（2）用MUSIC算法，捕捉图像的两个尖峰，得到频率的估计值：
图2.1MUSIC算法得到的谱峰
多运行几次，得到的频率值如下：
，
，
，
可以看出，MUSIC算法基本能够分辨出0.5和0.52两个频率。
相比于PHD算法，其误差小一些，估计的更准一些。
3附录
3.1题目1源程序
%%2阶自回归模型参数估计LSM&LSL算法by刘航20151108
clear;clf;clc;
%%初始化参数
%N是信号点数；
N=500;
na=1:
N;
%二阶自回归信号模型
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）;
fori=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;gama（1:
M,1）=1.0;
sef（1:
M,1）=1;seb（1:
M,1）=1;%前向和后向预测误差初始值
form=1:
M-1%按阶数迭代（阶数为2,matlab下标从1开始，故M设为3）
forn=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）-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）=Deta（m+1,n）/seb（m,n-1）;
end
end
kb（:
1）=0;kf（:
1）=0;
forn=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）,holdon;
plot（na,w（2,na）,'b--','linewidth',2）,holdon;
plot（na,a1n,'g--','linewidth',2）;holdon;
plot（na,a2n,'r--','linewidth',2）;
legend（'表示LMS中a1','表示LMS中a2','表示LSL中a1','表示LSL中a2',0）;holdon;
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性能参数比较'）;
gridon;
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）,holdon;
plot（na,a1_1n,'r--','linewidth',2）,holdon;
plot（na,a1_2n,'b--','linewidth',2）,holdon;
legend（'1.0','10','0.1',0）;
gridon;
%}
3.2题目2源程序
%%PHD算法by刘航20151108
%%初始化参数，产生复白噪声噪声和待估计信号
clc;
clearall;
N_sample=25;%取样自相关点数
M=2;%正弦信号的个数
sigma=0.1;
v_real=sigma*randn（N_sample,1）/sqrt
（2）;%白噪声的实部
v_imag=sigma*randn（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_sample:
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）;%零点的角度
fori=1:
M
ifw（i）<0
w（i）=w（i）+2*pi%angle函数求出的角度范围是[-pi,pi）,当求出的角频
%率w<０时，需要加上一个2*pi以调整到[0,2*pi）之间
end
end
f_PHD=w/（2*pi）
%%MUSCI算法by刘航20151109
%%初始化
clc;
clearall;
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）/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.52*n+pi/4））;
x=s+v;%观察信号
%%数据自相关矩阵获得；
x_corr=xcorr（x,'unbiased'）;%相关序列的无偏估计unbiased
Rx=（toeplitz（x_corr（N_sample:
N_sample+N-1）））.';%求自相关矩阵
%%求Rx的特征值和特征向量
[V,D]=eig（Rx）;%求特征值和特征向量
vmin=V（:
1:
N-M）;%噪声子空间的N-M个特征向量
%%通过方程（4.155）用作图法求估计频率值
K=10000;%在0--1之间分10000步搜索
f=linspace（0,1,K）;
e=zeros（N,K）;
forn=1:
N
e（n,:
）=exp（1i*（n-1）*2*pi*f）;
end
y=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
升级会员 