ADSP仿真报告 华中科技大学Word格式.docx
《ADSP仿真报告 华中科技大学Word格式.docx》由会员分享,可在线阅读,更多相关《ADSP仿真报告 华中科技大学Word格式.docx(21页珍藏版)》请在冰点文库上搜索。
为克服以上缺点,引入“最小二乘(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:
%前向和后向预测误差初始值
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
kb(:
1)=0;
kf(:
forn=1:
a1n(n)=kb(2,n)-kf(2,n)*kb(3,n);
a2n(n)=kb(3,n);
%%画图程序
figure
(1),plot(na,w(1,na),'
k--'
'
linewidth'
2),holdon;
plot(na,w(2,na),'
b--'
plot(na,a1n,'
g--'
2);
holdon;
plot(na,a2n,'
r--'
legend('
表示LMS中a1'
表示LMS中a2'
表示LSL中a1'
表示LSL中a2'
0);
plot(1.558*ones(n),'
y--'
);
plot(-0.81*ones(n),'
text(0,1.558,'
1.558'
color'
r'
),text(0,-0.81,'
-0.81'
title('
LSL与LMS性能参数比较'
);
gridon;
xlabel('
n'
ylabel('
a1(n),a2(n)'
%%LSL不同初始值收敛性能
%{
%更改sef(1:
M,1)=1的值为0.1和10(得到a1_1n,a1_2n)
figure
(2),plot(na,a1n,'
plot(na,a1_1n,'
plot(na,a1_2n,'
1.0'
10'
0.1'
%}
3.2题目2源程序
%%PHD算法by刘航20151108
%%初始化参数,产生复白噪声噪声和待估计信号
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<
0时,需要加上一个2*pi以调整到[0,2*pi)之间
f_PHD=w/(2*pi)
%%MUSCI算法by刘航20151109
%%初始化
N=12;
%题(5)中自相关矩阵Rx的阶数
%%产生复白噪声噪声和信号
%白噪声功率
SNR=10*log10(1/(sigma.^2));
%求信噪比
s=exp(1i*2*pi*0.5*n)+exp(1i*(2*pi*0.52*n+pi/4));
%观察信号
%%数据自相关矩阵获得;
x_corr=xcorr(x,'
%相关序列的无偏估计unbiased
%求自相关矩阵
%%求Rx的特征值和特征向量
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);
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