《现代数字信号处理及其应用》仿真作业.docx
- 文档编号:12220841
- 上传时间:2023-04-17
- 格式:DOCX
- 页数:59
- 大小:666.86KB
《现代数字信号处理及其应用》仿真作业.docx
《《现代数字信号处理及其应用》仿真作业.docx》由会员分享,可在线阅读,更多相关《《现代数字信号处理及其应用》仿真作业.docx(59页珍藏版)》请在冰豆网上搜索。
《现代数字信号处理及其应用》仿真作业
第三章仿真作业
3.17
(1)
N=32;
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(1i*2*pi*f1*(0:
N-1));
signal2=A2*exp(1i*2*pi*f2*(0:
N-1));
signal3=A3*exp(1i*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');
r11=real(r1);
r12=imag(r1);
r1=real(r);
r2=imag(r);
m=1-N:
N-1;
subplot(2,2,1);
stem(m,r11,'o');
xlabel('m');
ylabel('实部');
title('基于FFT的自相关函数快速计算');
subplot(2,2,2);
stem(m,r12,'o');
xlabel('m');
ylabel('虚部');
subplot(2,2,3);
stem(m,r1);
xlabel('m');
ylabel('实部');
title('教材式(3.1.2)估计的自相关函数');
subplot(2,2,4);
stem(m,r2);
xlabel('m');
ylabel('虚部');
(2)
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(1i*2*pi*f1*(0:
N-1));
signal2=A2*exp(1i*2*pi*f2*(0:
N-1));
signal3=A3*exp(1i*2*pi*f3*(0:
N-1));
un=signal1+signal2+signal3+noise;
NF=1024;
Spr=fftshift((1/NF)*abs(fft(un,NF)).^2);
Sprmax=max(Spr);
Spr=Spr/Sprmax;
f=(-(NF/2)+1:
(NF/2))/NF;
plot(f,20*log(Spr));
xlabel('f');
ylabel('归一化功率谱/dB');
title('周期图法');
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(1i*2*pi*f1*(0:
N-1));
signal2=A2*exp(1i*2*pi*f2*(0:
N-1));
signal3=A3*exp(1i*2*pi*f3*(0:
N-1));
un=signal1+signal2+signal3+noise;
M=64;
r=xcorr(un,N-1,'biased');
NF=1024;
BT=fftshift(fft(r,NF));
BTmax=max(BT);
BT=BT/BTmax;%归一化幅度
f=(-(NF/2)+1:
(NF/2))/NF;
plot(f,20*log(BT));
xlabel('f');
ylabel('归一化功率谱/dB');
title('BT法');
(3)
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(1i*2*pi*f1*(0:
N-1));
signal2=A2*exp(1i*2*pi*f2*(0:
N-1));
signal3=A3*exp(1i*2*pi*f3*(0:
N-1));
un=signal1+signal2+signal3+noise;
p=16;
r0=xcorr(un,p,'biased');
r=r0(p+1:
2*p+1);%从p+1开始取
a(1,1)=-r
(2)/r
(1);
sigma
(1)=r
(1)-(abs(r
(2)).^2)/r
(1);
form=2:
p
k(m)=-(r(m+1)+sum(a(m-1,1:
m-1).*r(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).^2));
end
NF=1024;
Par=sigma(p)./fftshift(abs(fft([1,a(p,:
)],NF)).^2);
Parmax=max(Par);
Par=Par/Parmax;%归一化幅度
f=(-(NF/2)+1:
(NF/2))/NF;
plot(f,20*log(Par));
xlabel('f');
ylabel('归一化功率谱/dB');
title('16阶AR模型的功率谱估计');
3.20
(1)
clc;
clearall;
closeall;
N=1000;
noise=(randn(1,N)+1i*randn(1,N))/sqrt
(2);
signal1=exp(j*0.5*pi*(0:
N-1)+j*unifrnd(0,2*pi,1,1));
signal2=exp(-j*0.3*pi*(0:
N-1)+j*unifrnd(0,2*pi,1,1));
un=signal1+signal2+noise;
M=8;
xs=zeros(M,N-M);
fork=1:
N-M
xs(:
k)=un(k+M-1:
-1:
k).';
end
R=(xs*xs')/(N-M);
[U,E,V]=svd(R);
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:
1,m);
end
z=roots(co);
ph=angle(z)/(2*pi);
z1=abs(z);
z2=abs(z1-1);
[estallv,estain]=sort(z2,'ascend');
所以单Root-MUSIC算法中最近单位圆的两个根为
0.0033-0.9977i0.5856+0.8074i
上述根的相位对应的归一化频率为
0.2495-0.1501
(2)
clc;clearall;closeall;
N=1000;
noise=(randn(1,N)+1i*randn(1,N))/sqrt
(2);
signal1=exp(j*0.5*pi*(0:
N-1)+j*2*pi*rand);
signal2=exp(-j*0.3*pi*(0:
N-1)+j*2*pi*rand);
un=signal1+signal2+noise;
M=8;
xs=zeros(M,N-M);
fork=1:
N-M
xs(:
k)=un(k+M-1:
-1:
k).';
end
R=(xs*xs')/(N-M);
[U,E,V]=svd(R);
G=U(:
3:
M);
f=[-0.5:
1/999:
0.5];
forff=1:
length(f)
w=f(ff)*2*pi;
forl=1:
M
aw(l)=exp(-j*w*(l-1));%计算a(w)
end
WW=aw*G*G'*aw';
Pmusic(ff)=abs(1./WW);%谱扫描函数
end
Pmusic=10*log10(Pmusic);
f=[-0.5:
1/999:
0.5];
plot(f,Pmusic);
holdon%绘图输出
xlabel('ω/2π')
ylabel('归一化Music谱/dB');
3.21
clc;
clearall;
closeall;
N=1000;
noise=(randn(1,N)+1i*randn(1,N))/sqrt
(2);
signal1=exp(j*0.5*pi*(0:
N-1)+j*unifrnd(0,2*pi,1,1));
signal2=exp(-j*0.3*pi*(0:
N-1)+j*unifrnd(0,2*pi,1,1));
un=signal1+signal2+noise;
M=8;
xs=zeros(M,N-M);
fork=1:
N-M
xs(:
k)=un(k+M-1:
-1:
k).';
end
Rxx=xs(:
1:
length(xs)-1)*xs(:
1:
length(xs)-1)'/(N-M-1);
Rxy=xs(:
1:
length(xs)-1)*xs(:
2:
length(xs))'/(N-M-1);
[U,E]=svd(Rxx);
emin=min(ev);
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);
ph=angle(z)/(2*pi);
z1=abs(z);
z2=abs(z1-1);
[estallv,estain]=sort(z2,'ascend');
所以单次ESPRIT算法中最近单位圆的两个特征值为
0.5845+0.8117i0.0036-0.9994i
上述特征值的相位对应的归一化频率为
-0.15070.2494
第四章仿真作业
4.18
%
(1)产生N=512点的样本序列
data_len=512;%样本序列的长度
trials=100;%随机试验的次数
A=zeros(data_len,2);EA=zeros(data_len,1);
B=zeros(data_len,2);EB=zeros(data_len,1);
form=1:
trials
a1=-0.975;
a2=0.95;
sigma_v_2=0.0731;
v=sqrt(sigma_v_2)*randn(data_len,1,trials);%产生v(n)
u0=[00];
num=1;
den=[1a1a2];
Zi=filtic(num,den,u0);%滤波器的初始条件
u=filter(num,den,v,Zi);%产生样本序列u(n)
%
(2)用LMS滤波器来估计w1和w2
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);
%LMS迭代过程
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
A=A+conj(w1)';
EA=EA+e1.^2;
B=B+conj(w2)';
EB=EB+e2.^2;
end
%剩余均方误差和失调参数
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
sJmin=sum(Jmin)/trials;
e1_100trials_ave=sum(e1)/trials;
e2_100trials_ave=sum(e2)/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
figure
(1);
plot(A/trials);holdon;
plot(conj(w1)');
xlabel('迭代次数');ylabel('权向量');title('步长为0.05权向量收敛曲线');
figure
(2);
plot(B/trials);holdon;
plot(conj(w2)');
xlabel('迭代次数');ylabel('权向量');title('步长为0.005权向量收敛曲线');
figure(3);
plot(EA/trials,'*');holdon;
plot(EB/trials,'-');
xlabel('迭代次数');ylabel('均方误差');title('步长分别为0.05和0.005学习曲线');
仿真结果:
剩余均方误差和失调参数结果分别为:
M1=0.0502
M2=0.0048
第五章仿真作业
5.10
解:
(1)求解M=2时的相关矩阵。
根据u(n)表达式可知a1=-0.99,σ2=0.93627,构造Yule-Walker方程
由于u(n)为是信号,故根据自相关函数的共轭对称性,有r
(1)=r(-1),求解方程可得
故滤波器权系数个数为M=2时的相关系数为
(2)根据
(1)中同样的方法,可得滤波器权系数个数为M=3时的相关系数为
(3)对相关矩阵进行特征分解
R2=[47.048746.5783;46.578347.0487];
R3=[47.048746.578346.1125;46.578347.048746.5783;46.112546.578347.0487];
eig_value2=eig(R2)
eig_value3=eig(R3)
eig_spread2=max(eig_value2)/min(eig_value2)
eig_spread3=max(eig_value3)/min(eig_value3)
(4)根据LMS算法均方误差收敛特性,可知步长因子应满足
当M=2时,步长因子应满足
,M=3时,有
,所以题目中给的
不符合要求,故在试验中M=2时,令
。
clc;
clearall;
M=2;%抽头个数
data_len=10000;
N=500;
n=1:
data_len;
a1=-0.99;
sigma_v=0.936271;
v=sqrt(sigma_v)*randn(data_len,1,N);
u0=0;
num=1;
den=[1a1];
Zi=filtic(num,den,u0);
u=filter(num,den,v,Zi);
u_len1=0.0006;
w1=zeros(M,data_len);
e1=zeros(data_len,1);
d1=zeros(data_len,1);
ww1=zeros(2*N,data_len);
%%%%%%%%%%%%%求学习曲线
fori=1:
N
forn=M+1:
data_len-1
w1(:
n+1)=w1(:
n)+u_len1*u(n-1:
-1:
n-M,:
i)*conj(e1(n));
d1(n+1)=w1(:
n+1)'*u(n:
-1:
n-M+1,:
i);
e1(n+1)=u(n+1,:
i)-d1(n+1);
ee1(n)=e1(n)*e1(n);
end
we1(i,:
)=ee1;
end
Ewe1=mean(we1);
plot(Ewe1)
holdon
xlabel('迭代次数');
ylabel('均方误差');
title('步长为0.001的500次独立实验的学习曲线');
第六章仿真作业
6.13
M=8;%阶数
N=1000;
noise=(randn(1,N)+j*randn(1,N))/sqrt
(2);
f1=0.1;
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(1i*2*pi*f1*(0:
N-1));
signal2=A2*exp(1i*2*pi*f2*(0:
N-1));
signal3=A3*exp(1i*2*pi*f3*(0:
N-1));
un=signal1+signal2+signal3+noise;
A=zeros(M,N-M+1);
forn=1:
N-M+1
A(:
n)=un(M+n-1:
-1:
n);
end
[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);
fork=1:
P
form=1:
M
a(m,k)=exp(-j*omega(k)*(m-1));
end
end
Pmvdr=zeros(size(omega));
fork=1:
P
Pmvdr(k)=1/(a(:
k)'*invphi*a(:
k));
end
Pmvdr=abs(Pmvdr/max(Pmvdr));
Pmvdr=10*log10(Pmvdr);
f=(-(P/2)+1:
(P/2))/P;
plot(f,Pmvdr);
xlabel('f');
ylabel('归一化功率谱/dB');
title('M=8时MVDR谱');
6.15
clc;
clearall;
M=2;%抽头个数
data_len=1000;
N=500;
a1=0.99;
sigma_v=0.995;
W=zeros(data_len,2);
fori=1:
N
v=sqrt(sigma_v)*randn(data_len,1);
num=1;
den=[1a1];
u0=zeros(length(den)-1,1);
Zi=filtic(num,den,u0);
u=filter(num,den,v,Zi);
n0=1;
b=u(n0+1:
data_len);
L=length(b);
u1=[zeros(M-1,1).',u.'].';
A=zeros(M,L);
fork=1:
L
A(:
k)=u1(M-1+k:
-1:
k);
end
delta=0.004;
lambda=0.98;
w=zeros(M,L+1);
epsilon=zeros(L,1);
P1=eye(M)/delta;
fork=1:
L-1
PIn=P1*A(:
k+1);
denok=lambda+A(:
k+1).'*PIn;
kn=PIn./denok;
epsilon(k+1)=b(k+1)-w(:
k).'*A(:
k+1);
w(:
k+1)=w(:
k)+kn.*conj(epsilon(k+1));
P1=P1/lambda-kn*A(:
k+1).'*P1/lambda;
epsilon2(k+1)=epsilon(k+1)*epsilon(k+1);
end
MSE(i,:
)=epsilon2;%也可以求和在求平均
W=W+conj(w)';
end
figure
(1);
MSE=mean(MSE);
plot(MSE);
xlabel('迭代次数');
ylabel('MSE');
title('500次重复独立实验的学习曲线');
figure
(2);
plot(W/500);
holdon
plot(conj(w)');
xlabel('迭代次数');
ylabel('权值');
title('权向量收敛曲线');
第七章仿真作业
7.14
Clear;
closeall;
N=3000;
sigma_v=0.0332;
v=sqrt(sigma_v)*randn(1,N);
%根据给定的AR模型产生u(n)序列
a1=1.6;
a2=-1.46;
a3=0.616;
a4=-0.1525;
u(1:
4)=0;
fori=1:
(N-4)
u(i+4)=a1*u(i+
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 现代数字信号处理及其应用 现代 数字信号 处理 及其 应用 仿真 作业