现代数字信号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
升级会员 