现代数字信号matlab处理仿真题.docx
- 文档编号:13703877
- 上传时间:2023-06-16
- 格式:DOCX
- 页数:43
- 大小:198.02KB
现代数字信号matlab处理仿真题.docx
《现代数字信号matlab处理仿真题.docx》由会员分享,可在线阅读,更多相关《现代数字信号matlab处理仿真题.docx(43页珍藏版)》请在冰点文库上搜索。
现代数字信号matlab处理仿真题
3.17
(1)相关函数
仿真代码:
A1=getAk(SNR1);
A2=getAk(SNR2);
A3=getAk(SNR3);%求得信号的幅度;
noise=randn(1,N)+j*randn(1,N);%构建高斯白噪声;
s1=getSk(A1,f1,N);
s2=getSk(A2,f2,N);
s3=getSk(A3,f3,N);;%产生3个复正弦信号
vn=s1+s2+s3+noise;
vk=fft(vn,2*N);%对v(n)补N个零,然后做2N点FFT
swk=((abs(vk)).^2)/N;%计算功率谱估计S(ω)
r0=ifft(swk);%对S(k)做ifft得到
r=[r0(N+2:
2*N),r0(1:
N)];%根据教程3.1.8式可得
r1=xcorr(vn,N-1,'biased');%直接计算自相关函数
%%%%%%%%%%%%%%%%%%%%%%%%%取序列实部,虚部%%%%%%
real_r=real(r);
imag_r=imag(r);
real_r1=real(r1);
imag_r1=imag(r1);
subplot(2,2,1);
stem(real_r);
xlabel('基于FFT的自相关函数的实部');
ylabel('实部');
subplot(2,2,2);
stem(imag_r);
xlabel('基于FFT的自相关函数的虚部');
ylabel('虚部');
subplot(2,2,3);
stem(real_r1);
ylabel('实部');
xlabel('估计的自相关函数的实部');
subplot(2,2,4);
stem(imag_r1);
xlabel('估计的自相关函数的虚实部');
ylabel('虚部');
functionAK=getAk(SNR)%求得幅度
%%%%%%%%%%%%%%%%%%%由SNR=10log(A^2/2*σ^2)%%%%%%%
AK=((10^(SNR/10))*2)^0.5;
functionSk=getSk(Ak,fk,N)
Sk=Ak*exp(j*2*pi*fk*(0:
N-1));
仿真波形:
(2)BT法和周期法估计
仿真程序:
clearall;
clc;
%设定N值可以改变抽样信号的点数,设定M值可以设定加窗的大小,设定N3可以补零,确定实际求fft的点数。
N=256;%样本观察长度
M=64;%加窗长度
f1=0.15;
f2=0.17;
f3=0.26;%定义3个归一化正弦频率
SNR1=30;
SNR2=30;
SNR3=27;%定义三个正弦信号信噪比
A1=getAk(SNR1);
A2=getAk(SNR2);
A3=getAk(SNR3);%求得信号的幅度;
noise=randn(1,N)+j*randn(1,N);%构建高斯白噪声;
s1=getSk(A1,f1,N);
s2=getSk(A2,f2,N);
s3=getSk(A3,f3,N);%产生3个复正弦信号
un=s1+s2+s3+noise;
%%%%%%%%%%%%下面采用周期法估计的频谱%%%%%%%%%%%%%
N2=2*N;
u2=zeros(1,N2);
u2(1:
N)=un;%对信号序列进行补零
Uw=fft(u2);%求DFT变换
Sw1=zeros(1,N2);
forw=1:
N2
Sw1(w)=((abs(Uw(w)))^2)/N;%计算功率谱
end
r0=zeros(1,N2);
r0=ifft(Sw1);
r=zeros(1,N2);
r(1:
N)=r0(N+1:
N2);
r(N+1:
N2)=r0(1:
N);
M=64;%加窗处理
forn=1:
2*M
r1(n)=r((N2-2*M)/2+n);
end%加窗之后的序列为r1
N3=256;%实际进行fft变换的点数。
r2=zeros(1,N3);
r2(1:
2*M)=r1;%将加窗之后的r1进行了补零得到的是r2.
Sw2=fft(r2);
Sw2=log10(abs(Sw2));%取模取对数
temp=zeros(1,N3/2);
temp=Sw2(1:
N3/2);
Sw2(1:
N3/2)=Sw2((N3/2+1):
N3);
Sw2((N3/2+1):
N3)=temp;%将0-2*pi的序列折换成-pi到pi的序列。
d=1/(N3-1);
x=zeros(1,N3);
fori=1:
N3
x(i)=-0.5+(i-1)*d;
end
Sw1=log10(abs(Sw1));%取模取对数
temp=zeros(1,N2/2);
temp=Sw1(1:
N2/2);
Sw1(1:
N2/2)=Sw1((N2/2+1):
N2);
Sw1((N2/2+1):
N2)=temp;%将0-2*pi的序列折换成-pi到pi的序列。
l=1/(N2-1);
y=zeros(1,N2);
fork=1:
N2
y(k)=-0.5+(k-1)*l;
end
subplot(1,2,1);
plot(x,Sw2);
xlabel('BT法');
ylabel('归一化的功率谱/dB');
subplot(1,2,2);
plot(y,Sw1);
xlabel('周期法');
ylabel('归一化的功率谱/dB');
仿真波形:
由仿真结果可以看出:
BT法和周期法在相同的阶数下,周期法比BT法分辩率更高,但是波动比BT法更大。
(3)Levinson-Durbin迭代法
仿真代码:
clearall;
clc;
%设定N值可以改变抽样信号的点数,设定M值可以设定加窗的大小,设定N3可以补零,确定实际求fft的点数。
N=256;%样本观察长度
p=16;%设定AR模型参数
f1=0.15;
f2=0.17;
f3=0.26;%定义3个归一化正弦频率
SNR1=30;
SNR2=30;
SNR3=27;%定义三个正弦信号信噪比
A1=getAk(SNR1);
A2=getAk(SNR2);
A3=getAk(SNR3);%求得信号的幅度;
noise=randn(1,N)+j*randn(1,N);%构建高斯白噪声;
s1=getSk(A1,f1,N);
s2=getSk(A2,f2,N);
s3=getSk(A3,f3,N);%产生3个复正弦信号
un=s1+s2+s3+noise;
r0=xcorr(un,p,'biased');%直接计算自相关函数
rp=r0(p+1:
2*p+1);
%%%%%%%%Lervinson_Durbin算法%%%%%%%%
a(1,1)=-rp
(2)/rp
(1);
e
(1)=rp
(1)-abs(rp
(2))^2/rp
(1);%计算σ1
b=0;
form=2:
p
k(m)=-(rp(m+1)+sum(a(m-1,1:
m-1).*rp(m:
-1:
2)))/e(m-1);
a(m,m)=k(m);
ford=1:
m-1
a(m,d)=a(m-1,d)+k(m)*conj(a(m-1,m-d));
end
e(m)=e(m-1)*(1-abs(k(m))^2);
end
%%%%计算16阶AR功率谱%%%%%%%%%%%%%%%%%%%%%%
NF=1024;%%%FFT点数
Sw=e(p)./(abs(fft([1,a(p,:
)],NF)).^2);
%%%%%%此下作用相当于函数fftshift
temp=zeros(1,NF/2);
temp=Sw(1:
NF/2);
Sw(1:
NF/2)=Sw((NF/2+1):
NF);
Sw((NF/2+1):
NF)=temp;%将0-2*pi的序列折换成-pi到pi的序列。
SAR=zeros(1,NF);
forw=1:
NF
SAR(w)=log10((abs(Sw(w))^2));
end
d=1/(NF-1);
x=zeros(1,NF);
fori=1:
NF
x(i)=-0.5+(i-1)*d;
end
plot(x,SAR);
xlabel('w/2π');
ylabel('归一化的功率谱/dB');
仿真波形:
3.20
(1)Root-MUSIC
仿真代码:
clearall;
clc;
N=1000;%样本个数
f1=0.5;%信号频率
f2=-0.3;
s1=getSk(f1,N);%产生第一个信号
s2=getSk(f2,N);%产生第二个信号
noise=(randn(1,N)+j*randn(1,N))/sqrt
(2);%构建高斯白噪声;
un=s1+s2+noise;%产生带噪声的信号
K=2;%信号源个数
M=8;%自相关矩阵阶数
R=getR(un,N,M);%根据教材3-5-18求自相关矩阵
%%%%%%%root_MUSIC进行频率估计%%%%%%%%%%%%%
[U,E]=svd(R);%矩阵U是特征值向量组成的矩阵,E是对角矩阵,对角矩阵的排列是从大到小
e=diag(E);
G=U(:
3:
M);
GG=G*G';
co=zeros(2*M-1,1);%由教材可知是关于变量z的2(M-1)次方程,共有2(M-1)个根
form=1:
M
co(m:
m+M-1)=co(m:
m+M-1)+GG(M:
-1:
1,m);
end
z=roots(co);%求多项式的跟
erro=abs(abs(z)-1);%除掉大于1的增根(方程式在此条件下变形的)
ph=angle(z)/(2*pi);
%挑选出最接近单位圆的N个根
disp(z)
disp(erro)
disp(ph)
functionaw=getaw(w,M)
aw=zeros(M,1);
forj=1:
M
aw(j)=exp(-w*(j-1)*i);
end
functionR=getR(x,N,M)
L=N-M+1;
tempx=zeros(M,1);
R=zeros(M,M);
forn=M:
N
forj=1:
M
tempx(j)=x(n-(j-1));
end
R=R+tempx*tempx';
end
R=R/L;
functionSk=getSk(fk,N)
Sk=exp(j*2*pi*fk*(0:
N-1)+j*2*pi*rand);%产生0到2π上均分布的随机相位信号
仿真结果:
w1=-0.3004w2=0.4996
(2)MUSIC算法
仿真程序:
clearall;
clc;
N=1000;%样本个数
f1=0.5;%信号频率
f2=-0.3;
s1=getSk(f1,N);%产生第一个信号
s2=getSk(f2,N);%产生第二个信号
noise=(randn(1,N)+j*randn(1,N))/sqrt
(2);%构建高斯白噪声;
un=s1+s2+noise;%产生带噪声的信号
K=2;%信号源个数
M=8;%自相关矩阵阶数
R=getR(un,N,M);%根据教材3-5-18求自相关矩阵
NF=2048;%定义扫描点数
%%%%%%%MUSIC进行谱峰收索%%%%%%%%%%%%%
[U,E]=svd(R);%矩阵U是特征值向量组成的矩阵,E是对角矩阵,对角矩阵的排列是从大到小
e=diag(E);
G=U(:
3:
M);
d=1/(NF-1);%求画图用的横坐标的间隔。
h=zeros(1,NF);
fori=1:
NF
h(i)=-0.5+(i-1)*d;
end
y=zeros(1,NF);
forj=1:
NF
w=h(j)*2*pi;
aw=getaw(w,M);
y(j)=log10(abs(1/((aw')*G*(G')*aw)));
end
plot(h,y);
xlabel('w/2π');
ylabel('归一化的功率谱/dB');
title('MUSIC谱估计')
仿真波形:
由仿真结果可知:
Root-MUSIC算法与MUSIC算法相比,两者误差都很小。
3.21ESPRIT算法
仿真代码:
clearall;
clc;
formatlong
N=1000;%样本个数
f1=0.5;%信号频率
f2=-0.3;
s1=getSk(f1,N);%产生第一个信号
s2=getSk(f2,N);%产生第二个信号
noise=(randn(1,N)+j*randn(1,N))/sqrt
(2);%构建高斯白噪声;
un=s1+s2+noise;%产生带噪声的信号
K=2;%信号源个数
M=8;%自相关矩阵阶数
[X,R]=corrmtx(un,M);
Rxx=R(1:
M,1:
M);
Rxy=R(1:
M,2:
M+1);%计算Rxy
%%%%%%%EPRIT进行频率估计%%%%%%%%%%%%%
[U,E]=svd(Rxx);
e=diag(E);
emin=e(end);%获取最小特征值
%%%构造Z矩阵%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Z=zeros(M,M);
fori=1:
M-1
Z(i,i+1)=1;
end
Cxx=Rxx-emin*eye(M);
Cxy=Rxy-emin*Z;
[U,E]=eig(Cxx,Cxy);
z=diag(E);
ph=angle(z)/(2*pi);%求所有特征值的相应归一化特征值
erro=abs(abs(z)-1);%与单位圆之间的距离
disp(ph);
disp(erro);
functionR=getR(x,N,M)
L=N-M+1;
tempx=zeros(M,1);
R=zeros(M,M);
forn=M:
N
forj=1:
M
tempx(j)=x(n-(j-1));
end
R=R+tempx*tempx';
end
R=R/L;
functionSk=getSk(fk,N)
Sk=exp(j*2*pi*fk*(0:
N-1)+j*2*pi*rand);%产生0到2π上均分布的随机相位信号
仿真结果:
w1=0.496818689408838w2=-0.285563358651823
4.18
仿真代码:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%程序功能:
产生512点的样本函数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clearall;
clc;
N=512;%样本序列长度
M=100;%独立试验次数
sigma_v_2=0.0731;
vn=sqrt(sigma_v_2)*randn(N,1,M);%产生均值为零,方差为0.0731的加性噪声,产生100次512X1的噪声用于100次独立试验
u0=[00];
a=[1-0.9750.95];
un=filter(1,a,vn);%构建AR模型,其中vn为输入,un为输出,产生100组独立信号1x512
%display(un);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%程序功能:
用LMS算法来估计w1,w2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u1=0.05;%定义两个步长
u2=0.005;
w1=zeros(2,N,M);%初始化权向量.两个列向量对应要估计的w1,w2100个2xN权向量
w2=zeros(2,N,M);%分别针对u1=0.05,u2=0.005
form=1:
M;%100次独立实验
e1(:
:
m)=zeros(N,1);%初始化相关参数
e2(:
:
m)=zeros(N,1);
d1(:
:
m)=zeros(N,1);
d2(:
:
m)=zeros(N,1);
forn=3:
N-1%信号向量时刻迭代
w1(:
n+1,m)=w1(:
n,m)+u1*un(n-1:
-1:
n-2,:
m)*conj(e1(n,1,m));
w2(:
n+1,m)=w2(:
n,m)+u2*un(n-1:
-1:
n-2,:
m)*conj(e2(n,1,m));
d1(n+1,1,m)=w1(:
n+1,m)'*un(n:
-1:
n-1,:
m);
d2(n+1,1,m)=w2(:
n+1,m)'*un(n:
-1:
n-1,:
m);
e1(n+1,1,m)=un(n+1,:
m)-d1(n+1,1,m);
e2(n+1,1,m)=un(n+1,:
m)-d2(n+1,1,m);
end
end
w1_ave=zeros(2,N);
w2_ave=zeros(2,N);
e1_ave=zeros(N,1);
e2_ave=zeros(N,1);
form=1:
M
w1_ave=w1_ave+w1(:
:
m);
w2_ave=w2_ave+w2(:
:
m);
e1_ave=e1_ave+e1(:
:
m).^2;
e2_ave=e2_ave+e2(:
:
m).^2;
end
w1_ave=w1_ave/M;%100次独立实验权向量的均值
w2_ave=w2_ave/M;
e1_ave=e1_ave/M;%100次独立实验的学习曲线均值
e2_ave=e2_ave/M;
t=1:
N;
figure
(1)
plot(t,w1(1,:
100),t,w1(2,:
100),t,w1_ave(1,:
),t,w1_ave(2,:
))
xlabel('迭代次数');ylabel('权向量')
title('步长=0.05')
figure
(2)
plot(t,w2(1,:
100),t,w2(2,:
100),t,w2_ave(1,:
),t,w2_ave(2,:
))
xlabel('迭代次数');ylabel('权向量')
title('步长=0.005')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%程序功能:
独立实验100次计算剩余误差和失调参数,画出学习曲线
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wopt=zeros(2,M);
Jmin=zeros(1,M);
sum_eig=zeros(M,1);
%%%%通过维纳-霍夫方程计算最小均方误差
%直接计算自相关函数
form=1:
M
rm=xcorr(un(:
:
m),'biased');
R=[rm(512),rm(513);rm(511),rm(512)];
p=[rm(512);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)/M;%100次平均误差
Jex1=e1_ave-sJmin;
Jex2=e2_ave-sJmin;%计算剩余均方误差
sum_eig_100=sum(sum_eig)/M;%100次特征值总和
Jexfin1=u1*sJmin*(sum_eig_100/(2-u1*sum_eig_100));
Jexfin2=u2*sJmin*(sum_eig_100/(2-u2*sum_eig_100));
M1=Jexfin1/sJmin;%失调参数m1
M2=Jexfin2/sJmin;%失调参数m2
figure(3)
plot(t,e1_ave,'*',t,e2_ave,'r');
xlabel('迭代次数');
ylabel('均方误差');
legend('u1=0.05','u2=0.005')
title('学习曲线');
axis([0,600,0,1])
display(M1);
display(M2);
仿真波形:
M1=-0.002061056056484
M2=-2.064886318291408e-04
由仿真结果可知:
当步长较大时,收敛速度较快,但失调和剩余均方误差较大,从而使稳态性能变差;而步长较小时,收敛速度虽然较慢,但失调和剩余均方误差减小,从而可以改善稳态性能;
问题:
当计算失调参数时,采用教材式(4.4.28)时,仿真得到的失调参数,较式(4.4.27)误差较大。
分析原因,式(4.4.28)只有在步长非常小时,误差才小。
%M1=u1*sum_eig_100/(2-u1*sum_eig_100);%由教材式(4.4.28)
%M2=u2*sum_eig_100/(2-u2*sum_eig_100);
5.10仿真结果及图形
(1)M=2时,
,求解Yule-Walker方程
得到相关矩阵
计算程序为:
%由yuler-walker方程计算M=2时的自相关矩阵
r2=inv([1,-0.99;-0.99,1])*[v_sigmal;0];
R2=[r2
(1),r2
(2);r2
(2),r2
(1)];%M1=2
(2)M=3时,
,求解Yule-Walker方程
得到自相关矩阵为
计算程序为:
r3=inv([1,-0.99,0;-0.99,1,0;0,-0.99,1])*[0.93627;0;0];
R3=[r3
(1),r3
(2),r3(3);r3
(2),r3
(1),r3
(2);r3(3),r3
(2),r3
(1)];%M2=3
(3)计算特征值扩展
%%%%%%%计算特征扩展值%%%%%%%%
eig1=eig(R2);%计算特征值
eig2=eig(R3);
eig_spread1=max(eig1)/min(eig1);
eig_spread2=max(eig2)/min(eig2);
%display(eig_spread1);
%display(eig_spread2);
M=2时特征值扩展是199.0000;
M=3时特征值扩展是444.2790。
(4)根据LMS算法均方误差收敛特性,M=2时步长因子应在区间(0,0.0106)之间,M=3时,步长因子应在区间(0,0.00708)之间,因此题中的步长因子不合理,仿真时不收敛。
故在仿真中,M=2时采用步长因子0.001,M=3时采用步长因子0.0006.
仿真程序:
%%%%%%%%%用LMS
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 现代 数字信号 matlab 处理 仿真