现代数字信号处理仿真作业.docx
《现代数字信号处理仿真作业.docx》由会员分享，可在线阅读，更多相关《现代数字信号处理仿真作业.docx（32页珍藏版）》请在冰豆网上搜索。
现代数字信号处理仿真作业
第三章仿真作业
3.17
（1）代码
clear;
N=32;
m=[-N+1:
N-1];
noise=（randn（1,N）+j*randn（1,N））/sqrt
（2）;
f1=0.15;
f2=0.17;
f3=0.26;
SNR1=30;
SNR2=30;
SNR3=27;
A1=10^（SNR1/20）;
A2=10^（SNR2/20）;
A3=10^（SNR3/20）;
signal1=A1*exp（j*2*pi*f1*（0:
N-1））;
signal2=A2*exp（j*2*pi*f2*（0:
N-1））;
signal3=A3*exp（j*2*pi*f3*（0:
N-1））;
un=signal1+signal2+signal3+noise;
uk=fft（un,2*N）;
sk=（1/N）*abs（uk）.^2;
r0=ifft（sk）;
r1=[r0（N+2:
2*N）,r0（1:
N）];
r=xcorr（un,N-1,'biased'）;
figure
subplot（2,2,1）
stem（m,real（r1））;
xlabel（'m'）;
ylabel（'FFT估计r1实部'）;
subplot（2,2,2）
stem（m,imag（r1））;
xlabel（'m'）;
ylabel（'FFT估计r1虚部'）;
subplot（2,2,3）
stem（m,real（r））;
xlabel（'m'）;
ylabel（'平均估计r实部'）;
subplot（2,2,4）
stem（m,imag（r））;
xlabel（'m'）;
ylabel（'平均估计r虚部'）;
仿真结果
（2）代码
clear;
N=256;
noise=（randn（1,N）+j*randn（1,N））/sqrt
（2）;
f1=0.15;
f2=0.17;
f3=0.26;
SNR1=30;
SNR2=30;
SNR3=27;
A1=10^（SNR1/20）;
A2=10^（SNR2/20）;
A3=10^（SNR3/20）;
signal1=A1*exp（j*2*pi*f1*（0:
N-1））;
signal2=A2*exp（j*2*pi*f2*（0:
N-1））;
signal3=A3*exp（j*2*pi*f3*（0:
N-1））;
un=signal1+signal2+signal3+noise;
NF=1024;
spr=fftshift（（1/NF）*abs（fft（un,NF））.^2）;
f1=（0:
length（spr）-1）*（1/（length（spr）-1））-0.5;
M=64;
r=xcorr（un,M,'biased'）;
bt=fftshift（abs（fft（r,NF）））;
f2=（0:
length（bt）-1）*（1/（length（bt）-1））-0.5;
figure
subplot（1,2,1）
plot（f1,10*log10（spr/max（spr）））;
xlabel（'w/2pi'）;
仿真结果
（3）
代码
clear;
N=1000;
fai1=rand（1,1）*2*pi;
fai2=rand（1,1）*2*pi;
noise=（randn（1,N）+j*randn（1,N））/sqrt
（2）;
un=exp（j*0.5*pi*（0:
N-1）+j*fai1）+exp（-j*0.3*pi*（0:
N-1）+j*fai2）+noise;
p=8;
cx=xcorr（un,p,'biased'）;
rxx=cx（p+1:
2*p）';
R=toeplitz（rxx）;
[u,s]=eig（R）;
nw=128;
ww=[-128:
128]/128*pi;
e=exp（-j*ww'*[0:
p-1]）%k行m列
ev=e*u（:
1:
p-2）;
pw=1./real（diag（ev*ev'））;
plot（ww/（2*pi）,10*log10（pw）/max（pw））;
仿真结果
3.20
（1）代码
clear;
N=1000;
fai1=rand（1,1）*2*pi;
fai2=rand（1,1）*2*pi;
noise=（randn（1,N）+j*randn（1,N））/sqrt
（2）;
un=exp（j*0.5*pi*（0:
N-1）+j*fai1）+exp（-j*0.3*pi*（0:
N-1）+j*fai2）+noise;
p=8;
cx=xcorr（un,p,'biased'）;
rxx=cx（p+1:
2*p）';
R=toeplitz（rxx）;
[u,s]=eig（R）;
nw=128;
ww=[-128:
128]/128*pi;
e=exp（-j*ww'*[0:
p-1]）%k行m列
ev=e*u（:
1:
p-2）;
pw=1./real（diag（ev*ev'））;
plot（ww/（2*pi）,10*log10（pw）/max（pw））;仿真结果
距离单位圆最近的两个解为-0.2363-0.9717i和0.3747+0.9271i，对应的归一化频率为0.1889和-0.2880
（2）代码
clear;
N=1000;
fai1=rand（1,1）*2*pi;
fai2=rand（1,1）*2*pi;
noise=（randn（1,N）+j*randn（1,N））/sqrt
（2）;
un=exp（j*0.5*pi*（0:
N-1）+j*fai1）+exp（-j*0.3*pi*（0:
N-1）+j*fai2）+noise;
p=8;
cx=xcorr（un,p,'biased'）;
rxx=cx（p+1:
2*p）';
R=toeplitz（rxx）;
[u,s]=eig（R）;
nw=128;
ww=[-128:
128]/128*pi;
e=exp（-j*ww'*[0:
p-1]）%k行m列
ev=e*u（:
1:
p-2）;
pw=1./real（diag（ev*ev'））;
plot（ww/（2*pi）,10*log10（pw）/max（pw））;
仿真结果
3.21
代码
clear;
N=1000;
fai1=rand（1,1）*2*pi;
fai2=rand（1,1）*2*pi;
noise=（randn（1,N）+j*randn（1,N））/sqrt
（2）;
un=exp（j*0.5*pi*（0:
N-1）+j*fai1）+exp（-j*0.3*pi*（0:
N-1）+j*fai2）+noise;
p=8;
fork=1:
N-p
xs（:
k）=un（k+p-1:
-1:
k）';
end
rxx=xs（:
1:
end-1）*xs（:
1:
end-1）'/（N-p-1）;
rxy=xs（:
1:
end-1）*xs（:
2:
end）'/（N-p-1）;
[u,e]=svd（rxx）;
ev=diag（e）;
emin=ev（end）;
z=[zeros（p-1,1）,eye（p-1）;0,zeros（1,p-1）];
cxx=rxx-emin*eye（p）;
cxy=rxy-emin*z;
[U,E]=eig（cxx,cxy）;
Z=diag（E）;
ph=angle（Z）/（2*pi）;
err=abs（abs（Z）-1）;
仿真结果
最接近单位圆的两个解分别为0.5867+0.8097i和0.0349-0.9984i，对应的归一化频率为0.1502和-0.2444。
第四章仿真题作业
4-18
（1）产生随机序列
代码
clear;
data_len=512;
trials=100;
n=1:
data_len;
a1=-0.975;
a2=0.95;
sigma_v_2=0.0731;
vn=sqrt（sigma_v_2）*randn（data_len,1,trials）;
u0=[00];
num=1;
den=[1a1a2];
zi=filtic（num,den,u0）;
u=filter（num,den,vn,zi）;
（2）单次实验估计最优权向量
代码
clear;
data_len=512;
trials=100;
n=1:
data_len;
a1=-0.975;
a2=0.95;
sigma_v_2=0.0731;
vn=sqrt（sigma_v_2）*randn（data_len,1,trials）;
u0=[00];
num=1;
m=1;
den=[1a1a2];
zi=filtic（num,den,u0）;
u=filter（num,den,vn,zi）;
mu1=0.05;
mu2=0.005;
w1=zeros（2,data_len）;
w2=zeros（2,data_len）;
e1=zeros（data_len,1）;
e2=zeros（data_len,1）;
d1=zeros（data_len,1）;
d2=zeros（data_len,1）;
%form=1:
trials
forn=3:
data_len-1
w1（:
n+1）=w1（:
n）+mu1*u（n-1:
-1:
n-2,:
m）*conj（e1（n））;
w2（:
n+1）=w2（:
n）+mu2*u（n-1:
-1:
n-2,:
m）*conj（e2（n））;
d1（n+1）=w1（:
n+1）'*u（n:
-1:
n-1,:
m）;
d2（n+1）=w2（:
n+1）'*u（n:
-1:
n-1,:
m）;
e1（n+1）=u（n+1,:
m）-d1（n+1）;
e2（n+1）=u（n+1,:
m）-d2（n+1）;
end
figure
（1）
holdon;
plot（w1（1,:
））;
plot（w1（2,:
））;
plot（w2（1,:
）,’g’）;
plot（w2（2,:
）,’g’）
仿真结果
（3）100次实验计算剩余均方误差
代码
clear;
data_len=512;
trials=100;
n=1:
data_len;
a1=-0.975;
a2=0.95;
sigma_v_2=0.0731;
vn=sqrt（sigma_v_2）*randn（data_len,1,trials）;
u0=[00];
num=1;
%m=1;
den=[1a1a2];
zi=filtic（num,den,u0）;
u=filter（num,den,vn,zi）;
mu1=0.05;
mu2=0.005;
w1=zeros（2,data_len）;
w2=zeros（2,data_len）;
e1=zeros（data_len,1）;
e2=zeros（data_len,1）;
d1=zeros（data_len,1）;
d2=zeros（data_len,1）;
jn1=zeros（1,trials）;
jn2=zeros（1,trials）;
mse1=zeros（1,trials）;
mse2=zeros（1,trials）;
form=1:
trials
forn=3:
data_len-1
w1（:
n+1）=w1（:
n）+mu1*u（n-1:
-1:
n-2,:
m）*conj（e1（n））;
w2（:
n+1）=w2（:
n）+mu2*u（n-1:
-1:
n-2,:
m）*conj（e2（n））;
d1（n+1）=w1（:
n+1）'*u（n:
-1:
n-1,:
m）;
d2（n+1）=w2（:
n+1）'*u（n:
-1:
n-1,:
m）;
e1（n+1）=u（n+1,:
m）-d1（n+1）;
e2（n+1）=u（n+1,:
m）-d2（n+1）;
end
jn1（m）=e1（256）;
jn2（m）=e2（256）
mse1（m）=abs（sum（jn1（1:
m））/m）;
mse2（m）=abs（sum（jn2（1:
m））/m）;
end
figure
（1）
holdon;
plot（mse1,'g'）;
plot（mse2,'b'）;
wopt=zeros（2,trials）;
jmin=zeros（1,trials）;
sum_eig=zeros（trials,1）;
form=1:
trials
rm=xcorr（u（:
:
m）,'biased'）;
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（:
m）;
sum_eig（m）=d（1,1）+d（2,2）;
end
figure
（2）
plot（jmin,'r'）;
sjmin=sum（jmin）/trials;
e1_100trials_ave=sum（jn1）/trials;
e2_100trials_ave=sum（jn2）/trials;
jex1=e1_100trials_ave-sjmin;
jex2=e2_100trials_ave-sjmin;
sum_eig_100trials=sum（sum_eig）/100;
jexfin=mu1*sjmin*（sum_eig_100trials/（2-mu1*sum_eig_100trials））;
jexfin2=mu2*sjmin*（sum_eig_100trials/（2-mu2*sum_eig_100trials））;
m1=jexfin/sjmin;
m2=jexfin2/sjmin;
仿真结果
LMS算法学习曲线
剩余均方误差jex1=-0.0569,jex2=-0.0959
失调参数：
m1=0.0523,m2=0.0050
第五章仿真题
1、当M=2时构造yule_walk方程
=
得到
故自相关矩阵
2、当M=3时构造Yule_walk方程
=
得到
故自相关矩阵为
3、r1=[47.048746.5783
46.578347.0487]
eig_value=eig（r）;
eig_spred=max（eig_value）/min（eig_value）;
r2=[47.048746.578346.1125
46.578347.048746.5783
46.112546.578347.0487]
eig_value=eig（r2）;
eig_spred=max（eig_value）/min（eig_value）;
特征值扩展结果
=199.0370,
=444.4107
4、代码
clear;
data_len=10000;
trials=100;
n=1:
data_len;
a1=-0.99;
M=2;
wsum=zeros（data_len,M）;
esum=zeros（1,data_len）;
sigma_v_2=0.93627;
form=1:
trials
vn=sqrt（sigma_v_2）*randn（data_len,1）;
u
（1）=vn
（1）;
fork=2:
data_len
u（k）=-a1*u（k-1）+vn（k）
end
w（1,:
）=zeros（1:
M）;
e
（1）=u
（1）;
mu=0.001;
uu=zeros（1,M）;
w（2,:
）=w（1,:
）+mu*e
（1）*uu;
uu=[u
（1）uu（1:
M-1）];
dd=（w（2,:
）*uu'）';
e
（2）=u（3）-dd;
fork=3:
data_len
w（k,:
）=w（k-1,:
）+mu*e（k-1）*uu;
uu=[u（k-1）uu（1:
M-1）];
dd=（w（k,:
）*uu'）';
e（k）=u（k）-dd;
end
wsum=wsum+w;
esum=esum+abs（e）;
end
figure
（1）
holdon;
plot（wsum/100）;
plot（w）;
figure
（2）
plot（esum/100）;
仿真结果
第六章仿真作业
6-13代码
clear;
m=4;
n=1000;
f1=[0.10.250.27];
SNR=[303027];
sigma=1;
am=sqrt（sigma*10.^（SNR/10））;
t=linspace（0,1,n）;
pha=2*pi*rand（size（f1））;
vn=sqrt（sigma/2）*randn（size（t））+j*sqrt（sigma/2）*randn（size（t））;
un=vn;
fork=1:
length（f1）
s=am（k）*exp（j*2*pi*n*f1（k）.*t+j*pha（k））;
un=un+s;
end%产生观测信号
A=zeros（m,n-m+1）;
fori=1:
n-m+1
A（:
i）=un（m+i-1:
-1:
i）;
end
[U,S,V]=svd（A'）;
invpha=V*inv（S'*S）*V';
p=1024;
f2=linspace（-0.5,0.5,p）;
omega=2*pi*f2;
a=zeros（m,p）;
fork=1:
p
fori=1:
m
a（i,k）=exp（-j*omega（k）*（i-1））;
end
end
pmvdr=zeros（size（omega））;
fork=1:
p
pmvdr（k）=1/（a（:
k）'*invpha*a（:
k））;
end
pmvdr=abs（pmvdr/max（abs（pmvdr）））;
plot（omega/（2*pi）,10*log10（pmvdr））;
仿真结果
6-15
a1=0.99;
sigma=0.995;
N=1000;
trials=500;
n0=1;
M=2;
L=N-n0;
vn=sqrt（sigma）*randn（N,1,trials）;
nume=1;
deno=[1a1];
u0=zeros（length（deno）-1,1）;
w=zeros（M,L+1,trials）;
epsilon=zeros（L,1,trials）;
W_mean=zeros（M,L+1）;
MSE=zeros（L,1）;
form=1:
trials
xic=filtic（nume,deno,u0）;
un=filter（nume,deno,vn（:
1,m）,xic）;
b=un（n0+1:
N）;
un1=[zeros（M-1,1）.',un.'].';
A=zeros（M,L）;
fork=1:
L
A（:
k）=un1（M-1+k:
-1:
k）;
end
delta=0.004;
lambda=0.98;
P1=eye（M）/delta;
fork=1:
L
PIn=P1*A（:
k）;
denok=lambda+A（:
k）'*PIn;
kn=PIn/denok;
epsilon（k,1,m）=b（k）-w（:
k,m）'*A（:
k）;
w（:
k+1,m）=w（:
k,m）+kn*conj（epsilon（k,1,m））;
P1=P1/lambda-kn*A（:
k）'*P1/lambda;
end
MSE=MSE+abs（epsilon（:
1,m））.^2;
W_mean=W_mean+w（:
:
m）;
end
figure
（1）;
plot（1:
L,1/trials*MSE）;
xlabel（'迭代次数'）;
ylabel（'MSE'）;
figure
（2）;
plot（1:
L+1,1/trials*W_mean）;
xlabel（'迭代次数'）;
ylabel（'权值'）;
第七章仿真作业
7-14代码
clearall;
N=3000;
gv=0.0332;
trials=100;
v=randn（N,1,trials）*sqrt（gv）;
a1=-2;
a2=1.4;
a3=-0.4;
a4=0.0384;
u1=zeros（N,1,trials）;
form=1:
trials
fori=1:
N-4
u1（i+4,1,m）=-a1*u1（i+3,1,m）-a2*u1（i+2,1,m）-a3*u1（i+1,1,m）...
-a4*u1（i,1,m）+v（i+4,1,m）;
end
end
N2=2000;
Jmin=0.005;
W_esti=zeros（4,N2,trials）;
W=zeros（4,N2+1）;
e=zeros（N2,1）;
form=1:
trials
P_esti=eye（4）;
forn=5:
N2
P_pre=P_esti;
U（:
n,m）=[u1（n-1,1,m）;u1（n-2,1,m）;u1（n-3,1,m）;u1（n-4,1,m）];
A=（U（:
n,m））'*P_pre*U（:
n,m）+Jmin;
K=P_pre*U（:
n,m）/A;
alpha（n,m）=u1（n,1,m）-（U（:
n,m））'*W_esti（:
n,m）;
W_esti（:
n+1,m）=W_esti（:
n,m）+K*alpha（n,m）;
P_esti=P_pre-K*（U（:
n,m））'*P_pre;
end
W=W+W_esti（:
:
m）;
e=e+abs（alpha（:
m））.^2;
end
W=1/trials*W;
e=1/trials*e;
h=1:
N2+1;
figure
（1）;
plot（h,W）;
xlabel（'迭代次数'）;
ylabel（'权值'）;
gtext（'w1'）;
gtext（'w2'）;
gtext（'w3'）;
gtext（'w4'）;
figure
（2）;
plot（1:
N2,e）;
xlabel（'迭代次数'）;
ylabel（'MSE'）;
第八章仿真题
clc;clear;
N=100;
M=10;
K=2;
theta=[-10;40]*pi/180;
SNR=[10;20];
Am=[sqrt（10.^（SNR/10））];
S=Am*ones（1,N）;
S（2,:
）=S（2,:
）.*exp（j*2*pi*rand（1,N））;
fora=1:
M
forb=1:
K
A（a,b）=exp（-j*（a-1）*pi*sin（theta（b）））;
end
end
V=zeros（M,N）;
form=1:
M
v=wgn（1,N,0,'complex'
升级会员 