现代数字信号处理仿真作业Word下载.docx
《现代数字信号处理仿真作业Word下载.docx》由会员分享，可在线阅读，更多相关《现代数字信号处理仿真作业Word下载.docx（26页珍藏版）》请在冰豆网上搜索。
r1=[r0（N+2:
2*N）,r0（1:
N）];
%%
r2=xcorr（un,N-1,'
biased'
）;
%画图
k=-N+1:
N-1;
figure
（1）
subplot（1,2,1）
stem（k,real（r1））
xlabel（'
m'
ylabel（'
实部'
subplot（1,2,2）
stem（k,imag（r1））
虚部'
figure
（2）
stem（k,real（r2））
stem（k,imag（r2））
%%周期图法
NF=1024;
Spr=fftshift（（1/NF）*abs（fft（un,NF））A2）;
kk=-0.5+（0:
NF-1）*（1/（NF-1））;
Spr_norm=10*log10（abs（Spr）/max（abs（Spr）））;
%%B祛
M=64;
r3=xcorr（un,M,'
BT=fftshift（fft（r3,NF））;
BT_norm=10*log10（abs（BT）/max（abs（BT）））;
figure（3）
plot（kk,Spr_norm）
w/2pi'
ylabel（title（'
周期图法'
）
subplot（1,2,2）plot（kk,BT_norm）
BT法'
）%%LDt代算法
归一化功率谱/DB'
p=16;
r0=xcorr（un,p,'
r4=r0（p+1:
2*p+1）;
%计算自相关函数
a（1,1）=-r4
（2）/r4
（1）;
sigma
（1）=r4
（1）-（abs（r4
（2）F2）/r4
（1）;
form=2:
p%LD^1代算法
k（m）=-（r4（m+1）+sum（a（m-1,1:
m-1）.*r4（m:
-1:
2）））/sigma（m-1）;
a（m,m）=k（m）;
fori=1:
m-1
a（m,i）=a（m-1,i）+k（m）*conj（a（m-1,m-i））;
end
sigma（m）=sigma（m-1）*（1-abs（k（m）F2）;
Par=sigma（p）./fftshift（abs（fft（[1,a（p,:
）],NF））.A2）;
Par_norm=10*log10（abs（Par）/max（abs（Par）））;
figure（4）
plot（kk,Par_norm）
归一化功率谱/DB'
title（'
16阶AR模型’）
%p阶AR模型的功率谱
2.仿真题3.20
单次Root-MUSIC算法中最接近单位圆的两个根为:
+
对应的归一化频率为：
相同信号的MUSIC谱估计结果如下
图对3.20信号进行MUSIC谱估计的结果
仿真程序（3_20）:
%%信号样本和高斯白噪声的产生
N=1000;
signal=[exp（1i*0.5*pi*（0:
N-1）+1i*2*pi*rand）;
%复正弦信号
exp（-1i*0.3*pi*（0:
N-1）+1i*2*pi*rand）];
%%计算自相关矩阵
M=8;
fork=1:
N-M
xs（:
k）=un（k+M-1:
k）.'
;
R=xs*xs'
/（N-M）;
%%自相关矩阵的特征值分解
[U,E]=svd（R）;
%%Root-MUSIC算法的实现
G=U（:
3:
M）;
Gr=G*G'
co=zeros（2*M-1,1）;
form=1:
M
co（m:
m+M-1）=co（m:
m+M-1）+Gr（M:
1,m）;
z=roots（co）;
ph=angle（z）/（2*pi）;
err=abs（abs（z）-1）;
%刻算MUSIC谱
En=U（:
2+1:
NF=2048;
forn=1:
NF
Aq=exp（-1i*2*pi*（-0.5+（n-1）/（NF-1））*（0:
M-1）'
）;
Pmusic（n）=1/（Aq'
*En*En'
*Aq）;
Pmusic_norm=10*log10（abs（Pmusic）/max（abs（Pmusic）））;
plot（kk,Pmusic_norm）
w/2*pi'
归一化功率谱/dB'
3.仿真题3.21
单次ESPRITT法中最接近单位元的两个特征值为：
-
仿真程序（3_21）:
%%自相关矩阵的计算
Rxx=xs（:
1:
end-1）*xs（:
end-1）'
/（N-M-1）;
Rxy=xs（:
2:
end）'
%%特征值分解
[U,E]=svd（Rxx）;
ev=diag（E）;
emin=ev（end）;
Z=[zeros（M-1,1）,eye（M-1）;
0,zeros（1,M-1）];
Cxx=Rxx-emin*eye（M）;
Cxy=Rxy-emin*Z;
%%广义特征值分解
[U,E]=eig（Cxx,Cxy）;
z=diag（E）;
4.仿真题4.18
步长为0.05时失调参数为m1=0.0493；
步长为0.005时失调参数为m2=0.0047。
图步长为0.05时权向量的收敛曲线
图步长为0.005时权向量的收敛曲线
图步长分别为0.05和0.005时100次独立实验的学习曲线
仿真程序（4_18）:
%%产生100组独立样本序列
data_len=512;
trials=100;
n=1:
data_len;
a1=-0.975;
a2=0.95;
sigma_v_2=0.0731;
v=sqrt（sigma_v_2）*randn（data_len,1,trials）;
u0=[00];
num=1;
den=[1,a1,a2];
Zi=filtic（num,den,u0）;
u=filter（num,den,v,Zi）;
%产生100组独立信号
%%LM链代
mu1=0.05;
mu2=0.005;
w1=zeros（2,data_len,trials）;
w2=w1;
100;
temp=zeros（data_len,1）;
e1（:
:
m）=temp;
e2（:
d1（:
d2（:
forn=3:
data_len-1
w1（:
n+1,m）=w1（:
n,m）+mu1*u（n-1:
n-2,:
m）*conj（e1（n,1,m））;
w2（:
n+1,m）=w2（:
n,m）+mu2*u（n-1:
m）*conj（e2（n,1,m））;
d1（n+1,1,m）=w1（:
n+1,m）'
*u（n:
n-1,:
m）;
d2（n+1,1,m）=w2（:
e1（n+1,1,m）=u（n+1,:
m）-d1（n+1,1,m）;
e2（n+1,1,m）=u（n+1,:
m）-d2（n+1,1,m）;
t=1:
w1_mean=zeros（2,data_len）;
w2_mean=w1_mean;
e1_mean=zeros（data_len,1）;
e2_mean=e1_mean;
100
w1_mean=w1_mean+w1（:
w2_mean=w2_mean+w2（:
e1_mean=e1_mean+e1（:
m）.A2;
e2_mean=e2_mean+e2（:
w1_mean=w1_mean/100;
%100次独立实验权向量的均值
w2_mean=w2_mean/100;
e1_100trials_ave=e1_mean/100;
%100次独立实验的学习曲线均值
e2_100trials_ave=e2_mean/100;
plot（t,w1（1,:
90）,t,w1（2,:
90）,t,w1_mean（1,:
）,t,w1_mean（2,:
））
迭代次数'
权向量'
步长=0.05'
plot（t,w2（1,:
90）,t,w2（2,:
90）,t,w2_mean（1,:
）,t,w2_mean（2,:
步长=0.005'
%%计算剩余误差和失调参数
wopt=zeros（2,trials）;
Jmin=zeros（1,trials）;
sum_eig=zeros（trials,1）;
trials
rm=xcorr（u（:
m）,'
R=[rm（512）,rm（513）;
rm（511）,rm（512）];
p=[rm（511）;
rm（510）];
wopt（:
m尸R\p;
[v,d]=eig（R）;
Jmin（m尸rm（512）-p'
*wopt（:
sum_eig（m）=d（1,1）+d（2,2）;
end
sJmin=sum（Jmin）/trials;
Jex1=e1_100trials_ave-sJmin;
%^余均方误差mu1
Jex2=e2_100trials_ave-sJmin;
%^余均方误差mu2
sum_eig_100trials=sum（sum_eig）/100;
Jexfin1=mu1*sJmin*（sum_eig_100trials/（2-mu1*sum_eig_100trials））;
Jexfin2=mu2*sJmin*（sum_eig_100trials/（2-mu2*sum_eig_100trials））;
M1=Jexfin1/sJmin;
%^调参数ml
M2=Jexfin2/sJmin;
%^调参数m2
plot（t,e1_100trials_ave,'
*'
t,e2_100trials_ave）
xlabel（'
迭代次数'
均方误差’）legend（'
u1=0.05'
'
u2=0.005'
axis（[0,600,0,1]）
5.仿真题5.10
（1）M=2时，a[=-0.99,a2=0.93627,求解Yule-Walker方程
可得到自相关矩阵相应的计算程序为r2=inv（[1,-0.99;
-0.99,1]）*[0.93627;
0];
R2=[r2
（1）,r2⑵;
r2
（2）,r2
（1）];
%M=2
（2）M=3时，a1=-0.99,a2=0户2=0.93627,求解Yule-Walker方程
可得到自相关矩阵为相应的计算程序为r3=inv（[1,-0.99,0;
-0.99,1,0;
0,-0.99,1]）*[0.93627;
0;
R3=[r3
（1）,r3
（2）,r3（3）;
r3
（2）,r3
（1）,r3
（2）;
r3（3）,r3
（2）,r3
（1）];
%M=3
（3）计算特征值扩展
%%寺征值分解
eig_value_1=eig（R2）;
eig_value_2=eig（R3）;
%%寺征值扩展
eig_spread_1=max（eig_value_1）/min（eig_value_1）;
eig_spread_2=max（eig_value_2）/min（eig_value_2）;
M=2时特征值扩展是199.0000;
M=3时特征值扩展是444.2790。
⑶根据LMS算法均方误差收敛特性，M=2时步长因子应在区间（0,0.0213）中，M=3时,步长因子应在区间（0,0.0142）之间，因此题中的步长因子不合理。
故在仿真中，M=2时
采用步长因子0.001,M=3时采用步长因子0.0006.
图500次独立实验M=2步长为0.001时权向量收敛曲线
图500次独立实验M=3步长为0.0006时权向量收敛曲线
图500次独立实验M=2步长为0.001时的学习曲线
图500次独立实验M=3步长为0.0006时的学习曲线
仿真程序（5_10_4）:
%%产生系统输入白噪声
L=10000;
sigma_v1_2=0.93627;
500
v（:
m）=sqrt（sigma_v1_2）*randn（L,1）;
%渔成500组独立的AR模型信号
a1=-0.99;
500u（1,1,m）=v（1,m）;
fork=2:
L
u（k,1,m）=-a1*u（k-1,1,m）+v（k,m）;
%%LM链代算法
M=2;
%M=3;
mu=0.001;
%mu=0.0006;
w=zeros（L,M,500）;
e（1,m）=u（1,m）;
uu=zeros（1,M）;
w（2,:
m）=w（1,:
m）+mu*e（1,m）*uu;
uu=[u（1,m）uu（1:
M-1）];
dd=（w（2,:
m）*uu'
）'
e（2,m）=u（3,m）-dd;
fork=3:
w（k,:
m）=w（k-1,:
m）+mu*e（k-1,m）*uu;
uu=[u（k-1,1,m）uu（1:
dd=（w（k,:
e（k,m）=u（k,m）-dd;
%%M=2
e_mean=zeros（10000,1）;
w_mean=zeros（10000,2）;
w_mean=w_mean+w（:
e_mean=e_mean+e（:
m『2;
w_mean=w_mean/500;
e_mean=e_mean/500;
L;
plot（t,w（:
1,100）,t,w（:
2,100）,t,w_mean（:
1）,t,w_mean（:
2））
迭代次数n'
抽头权值'
M=2,步长0.001的权向量收敛曲线'
figure
（2）plot（t,e_mean）
MSE'
M=2,步长0.001的学习曲线'
%%M=3
w_mean=zeros（10000,3）;
form=1:
2,100）,t,w（:
3,100）,t,w_mean（:
2）,t,w_m
ean（:
3））
M=3,步长0.0006的权向量收敛曲线'
M=2,步长0.0006的学习曲线'
6.仿真题6.13
滤波器抽头个数为4时
图M=4时的MVDR谱
图M=4时基于奇异值分解的MVDR谱
从上面两图可以看出，M=4时并没有将3个频点分辨出来，增加滤波器阶数可以解决此问
题，因此当M=20时仿真结果如下两图所示：
图M=20时的MVDR谱
图M=20时基于奇异值分解的MVDR谱
仿真程序（6_13）:
%%产生观测信号
M=4;
%M=20;
f=[0.10.250.27];
sigma=1;
Am=sqrt（sigma*10.A（SNR/10））;
t=linspace（0,1,N）;
phi=2*pi*rand（size（f））;
vn=sqrt（sigma/2）*randn（size（t））+1i*sqrt（sigma/2）*randn（size（t））;
Un=vn;
length（f）
s=Am（k）*exp（1i*2*pi*N*f（k）.*t+1i*phi（k））;
Un=Un+s;
Un=Un.'
%%构建矩阵
A=zeros（M,N-M+1）;
N-M+1
A（:
n）=Un（M+n-1:
n）;
[U,S,V]=svd（A'
invphi=V*inv（S'
*S）*V'
%%构建向量
P=1024;
f=linspace（-0.5,0.5,P）;
omega=2*pi*f;
a=zeros（M,P）;
P
a（m,k）=exp（-1i*omega（k）*（m-1））;
%刻算MVDRlf
Pmvdr=zeros（size（omega））;
Pmvdr（k）=1/（a（:
k）'
*invphi*a（:
k））;
Pmvdr=abs（Pmvdr/max（abs（Pmvdr）））;
Pmvdr=10*log10（Pmvdr）;
P-1）*（1/（P-1））;
plot（kk,Pmvdr）
%履于习题6.11的奇异值分解的MVD昉法
Sw=zeros（1,M）;
Sw（i）=（a（:
）*V（:
i）/S（i,i）;
P_svd（k）=1/sum（（abs（Sw））-2）;
P_svd=abs（P_svd/max（abs（P_svd）））;
P_svd=10*log10（P_svd）;
归一化频谱/dB'
M=4的MVDRif'
）
%title（'
M=20的MVDRf'
plot（kk,P_svd）
M=4的基于SVD的MVDR®
谱'
M=20的基于SVD的MVD砌谱'
7.仿真题6.15
图单次实验估计权值以及500次独立实验的估计权值
图500次独立实验的学习曲线
仿真程序（6_15）:
%花生AR模型的输入信号
a1=0.99;
sigma=0.995;
trials=500;
vn=sqrt（sigma）*randn（N,1,trials）;
nume=1;
deno=[1a1];
u0=zeros（length（deno）-1,1）;
xic=filtic（nume,deno,u0）;
un=filter（nume,deno,vn,xic）;
%产生500组独立的信号
%%产生期望信号和观测数据矩阵
n0=1;
b=un（n0+1:
N,:
L=length（b）;
A=zeros（M,L,trials）;
un1=[zeros（M-1,1）.'
un（:
m
升级会员 