现代数字信号处理仿真作业.docx
- 文档编号:94331
- 上传时间:2022-10-02
- 格式:DOCX
- 页数:31
- 大小:1.14MB
现代数字信号处理仿真作业.docx
《现代数字信号处理仿真作业.docx》由会员分享,可在线阅读,更多相关《现代数字信号处理仿真作业.docx(31页珍藏版)》请在冰豆网上搜索。
现代数字信号处理仿真作业
1.仿真题3.17
仿真结果及图形:
图1基于FFT的自相关函数计算
图2基于式3.1.2的自相关函数的计算
图3周期图法和BT法估计信号的功率谱
图4利用LD迭代对16阶AR模型的功率谱估计
16阶AR模型的系数为:
a1=-0.402637623107952-0.919787323662670i;
a2=-0.013530139693503+0.024214641171318i;
a3=-0.074241889634714-0.088834852915013i;
a4=0.027881022353997-0.040734794506749i;
a5=0.042128517350786+0.068932699075038i;
a6=-0.0042799971761507+0.028686095385146i;
a7=-0.048427890183189-0.019713457742372i;
a8=0.0028768633718672-0.047990801912420i
a9=0.023971346213842+0.046436389191530i;
a10=0.026025963987732+0.046882756497113i;
a11= -0.033929397784767-0.0053437929619510i;
a12=0.0082735406293574-0.016133618316269i;
a13=0.031893903622978-0.013709547028453i ;
a14=0.0099274520678052+0.022233240051564i;
a15=-0.0064643069578642+0.014130696335881i;
a16=-0.061704614407581-0.077423818476583i.
仿真程序(3_17):
clearall
clc
%%产生噪声序列
N=32;%基于FFT的样本长度
%N=256;%周期图法,BT法,AR模型功率谱估计的长度
vn=(randn(1,N)+1i*randn(1,N))/sqrt
(2);
%%产生复正弦信号
f=[0.150.170.26];%归一化频率
SNR=[303027];%信噪比
A=10.^(SNR./20);%幅度
signal=[A
(1)*exp(1i*2*pi*f
(1)*(0:
N-1));%复正弦信号
A
(2)*exp(1i*2*pi*f
(2)*(0:
N-1));
A(3)*exp(1i*2*pi*f(3)*(0:
N-1))];
%%产生观察样本
un=sum(signal)+vn;
%%利用3.1.1的FFT估计
Uk=fft(un,2*N);
Sk=(1/N)*abs(Uk).^2;
r0=ifft(Sk);
r1=[r0(N+2:
2*N),r0(1:
N)];
%%利用3.1.2估计R
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))
xlabel('m');ylabel('虚部');
figure
(2)
subplot(1,2,1)
stem(k,real(r2))
xlabel('m');ylabel('实部');
subplot(1,2,2)
stem(k,imag(r2))
xlabel('m');ylabel('虚部');
%%周期图法
NF=1024;
Spr=fftshift((1/NF)*abs(fft(un,NF)).^2);
kk=-0.5+(0:
NF-1)*(1/(NF-1));
Spr_norm=10*log10(abs(Spr)/max(abs(Spr)));
%%BT法
M=64;
r3=xcorr(un,M,'biased');
BT=fftshift(fft(r3,NF));
BT_norm=10*log10(abs(BT)/max(abs(BT)));
figure(3)
subplot(1,2,1)
plot(kk,Spr_norm)
xlabel('w/2pi');ylabel('归一化功率谱/DB');
title('周期图法')
subplot(1,2,2)
plot(kk,BT_norm)
xlabel('w/2pi');ylabel('归一化功率谱/DB');
title('BT法')
%%LD迭代算法
p=16;
r0=xcorr(un,p,'biased');
r4=r0(p+1:
2*p+1);%计算自相关函数
a(1,1)=-r4
(2)/r4
(1);
sigma
(1)=r4
(1)-(abs(r4
(2))^2)/r4
(1);
form=2:
p%LD迭代算法
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))^2);
end
Par=sigma(p)./fftshift(abs(fft([1,a(p,:
)],NF)).^2);%p阶AR模型的功率谱
Par_norm=10*log10(abs(Par)/max(abs(Par)));
figure(4)
plot(kk,Par_norm)
xlabel('w/2pi');ylabel('归一化功率谱/DB');
title('16阶AR模型')
2.仿真题3.20
仿真结果及图形:
单次Root-MUSIC算法中最接近单位圆的两个根为:
-0.001156047541561+1.001503153449793i
0.587376604261220-0.810845628739986i
对应的归一化频率为:
0.250183714447964
-0.150223406926494
相同信号的MUSIC谱估计结果如下
图5对3.20信号进行MUSIC谱估计的结果
仿真程序(3_20):
clearall
clc
%%信号样本和高斯白噪声的产生
N=1000;
vn=(randn(1,N)+1i*randn(1,N))/sqrt
(2);
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)];
un=sum(signal)+vn;
%%计算自相关矩阵
M=8;
fork=1:
N-M
xs(:
k)=un(k+M-1:
-1:
k).';
end
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:
1,m);
end
z=roots(co);
ph=angle(z)/(2*pi);
err=abs(abs(z)-1);
%%计算MUSIC谱
En=U(:
2+1:
M);
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);
end
kk=-0.5+(0:
NF-1)*(1/(NF-1));
Pmusic_norm=10*log10(abs(Pmusic)/max(abs(Pmusic)));
plot(kk,Pmusic_norm)
xlabel('w/2*pi');ylabel('归一化功率谱/dB')
3.仿真题3.21
仿真结果及图形:
单次ESPRIT算法中最接近单位元的两个特征值为:
0.001826505974929+1.000690248438859i
0.586994191014025-0.809491260856630i
对应的归一化频率为:
0.249709503383161
-0.150146235268272
仿真程序(3_21):
clearall
clc
%%信号样本和高斯白噪声的产生
N=1000;
vn=(randn(1,N)+1i*randn(1,N))/sqrt
(2);
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)];
un=sum(signal)+vn;
%%自相关矩阵的计算
M=8;
fork=1:
N-M
xs(:
k)=un(k+M-1:
-1:
k).';
end
Rxx=xs(:
1:
end-1)*xs(:
1:
end-1)'/(N-M-1);
Rxy=xs(:
1:
end-1)*xs(:
2:
end)'/(N-M-1);
%%特征值分解
[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);
ph=angle(z)/(2*pi);
err=abs(abs(z)-1);
4.仿真题4.18
仿真结果及图形:
步长为0.05时失调参数为m1=0.0493;
步长为0.005时失调参数为m2=0.0047。
图6步长为0.05时权向量的收敛曲线
图7步长为0.005时权向量的收敛曲线
图8步长分别为0.05和0.005时100次独立实验的学习曲线
仿真程序(4_18):
clearall
clc
%%产生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组独
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 现代 数字信号 处理 仿真 作业