信号处理仿真题作业.docx
《信号处理仿真题作业.docx》由会员分享，可在线阅读，更多相关《信号处理仿真题作业.docx（15页珍藏版）》请在冰点文库上搜索。
信号处理仿真题作业
信
号
处
理
基
础
仿
真
作
业
学号：
S1*******7
姓名：
贾雪婷
3.17在计算机上用如下方法产生随机信号
的观测样本：
首先产生一段零均值、方差为
的复高斯白噪声序列
；然后在
上叠加三个复正弦信号，它们的归一化频率分别是f1=0.15，f2=0.17和f3=0.26。
调整
和正弦信号的幅度，使在f1、f2和f3处得信噪比分别为30dB、30dB和27dB。
（1）令信号观测样本长度N=32，试用3.1.1节讨论的基于FFT的自相关函数快速计算方法估计出自相关函数
，并与教材式（3.1.2）估计出的自相关函数
做比较。
产生零均值、方差为1的复高斯白噪声序列y
>>y=randn（1,32）;
>>y=y-mean（y）;
>>y=y/std（y）;
>>a=0;
>>b=sqrt
（2）;
>>y=a+b*y
产生三个复正弦信号并产生观察样本：
>>N=32;
>>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+y
基于FFT的自相关函数快速计算方法：
N=32;
>>Uk=fft（un,2*N）;
Sk=（1/N）*abs（Uk）.^2;
r0=ifft（Sk）;
r1=[r0（N+2:
2*N）,r0（1:
N）];
>>figure
（1）;
>>stem（real（r1））;
>>figure
（2）;
>>stem（imag（r1））
输出结果为：
图1基于FFT的自相关函数快速计算
实部：
虚部：
教材中式（3.1.2）估计自相关函数
>>r=xcorr（un,N-1,'biased'）;
>>figure
（1）;
>>stem（real（r））
>>figure
（2）;
>>stem（imag（r））
输出结果为：
图2教材式（3.1.2）估计的自相关函数
实部：
虚部：
（2）令信号观测样本长度N=256，试用BT法和周期图法估计
的功率谱，这里设BT法中所用自相关函数的单边长度M=64。
周期图法
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））;
>>xn=signal1+signal2+signal3+noise
NF=1024;
Spr=fftshift（（1/NF）*abs（fft（xn,NF））.^2）;
Spr=Spr/max（Spr）;
t=[-0.5:
1/NF:
0.5-1/NF];
plot（t,10*log10（Spr））;
图3周期图法所得功率谱
BT法
>>M=64;
>>r=xcorr（xn,M,'biased'）;
>>NF=1024;
>>BT=fftshift（fft（r,NF））;
BT=BT/max（BT）;
>>plot（（-511:
512）/1024,10*log10（BT））图4BT法所得功率谱
（3）令信号观测样本长度N=256，试用Levinson-Durbin迭代算法求解AR模型的系数并估计
的功率谱，模型的阶数取为p=16。
p=16;
r0=xcorr（xn,p,'biased'）;
r=r0（p+1:
2*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）;
Par=Par/max（Par）;
plot（（-511:
512）/1024,10*log10（Par））
图516阶AR模型的功率谱估计
3.18设随机过程
为
其中，
是零均值、方差为1的白噪声，
、
是相互独立并在
上服从均匀分布的随机相位。
使用MVDR方法进行信号频率估计的方针实验，画出频率估计谱线，并给出正弦信号频率的估计值。
（要求：
信号样本数取1000，估计的自相关矩阵为8阶）
解：
clearall;
%产生0均值，方差为1的复高斯白噪声序列v（n）
N=1000;
FS=1;
noise=（randn（1,N）+j*randn（1,N））/sqrt
（2）;
%产生带噪声的信号样本u（n）
signal1=exp（j*0.5*pi*（0:
N-1）/FS+j*2*pi*rand）;
signal2=exp（-j*0.3*pi*（0:
N-1）/FS+j*2*pi*rand）;
un=signal1+signal2+noise;
%利用un来估计自相关函数
r=zeros（1,N）;
form=1:
N
forn=m:
N
r（m）=r（m）+un（n）*（un（n-（m-1）））';
end
r（m）=r（m）/N;
end
M=8;%设定横向滤波器的阶数
%利用自相关矩阵和自相关函数的关系，构建自相关矩阵
R=zeros（M,M）;
fori=1:
M
forj=1:
M
ifi<=j
R（i,j）=r（1+j-i）;%得到了M阶的矩阵R
else
R（i,j）=（r（1+i-j））';
end
end
end
N3=2048;%设定画图时描点的数目。
d=1/（N3-1）;%求画图用的横坐标的间隔。
h=zeros（1,N3）;
fori=1:
N3
h（i）=-0.5+（i-1）*d;
end
y=zeros（1,N3）;
forj=1:
N3
w=h（j）*2*pi;
y（j）=10*log10（abs（getPMVDR（w,M,R）））;
end
plot（h,y）;
title（'MVDR谱估计'）;
xlabel（'\omega/2\pi'）;ylabel（'归一化功率谱/dB'）;
axis（[-0.50.5-161]）;
图6MVDR谱估计
3.19复正弦加白噪声随机过程
同题3.18中所给。
试用MUSIC方法进行信号频率估计的仿真实验。
要求：
信号样本数取1000，估计的自相关矩阵为8阶）
（1）采用AIC和MDL准则估计信号源个数；
N=1000;
noise=（randn（1,N）+j*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;
fork=1:
N-M
xs（:
k）=un（k+M-1:
-1:
k）.';
end
R=xs*xs'/（N-M）;%自相关矩阵的特征值分解
[U,E]=svd（R）;%U是特征矢量组成的矩阵，E是对角阵，对角元素是由大到小排列的特征值
ev=diag（E）;%提取对角元素上的特征值；%根据AIC准则进行信号源个数的估计
fork=1:
M
dec=prod（ev（k:
M）.^（1/（M-k+1）））;%计算第一项中对数的自变量的分子
nec=mean（ev（k:
M））;%计算第一项中对数的自变量的分母
lnv=（dec/nec）^（（M-k+1）*N）;%计算第一项中对数的自变量
AIC（k）=-2*log（lnv）+2*（k-1）*（2*M-k+1）;%计算AIC准则
end
[Amin,K]=min（AIC）;
N1=K-1;%信号源个数估计
%根据MDL准则进行信号源个数的估计
fork=1:
M
dec=prod（ev（k:
M）.^（1/（M-k+1）））;
nec=mean（ev（k:
M））;
lnv=（dec/nec）^（（M-k+1）*N）;
MDL（k）=-log（lnv）+（k-1）/2*（2*M-k+1）*log（N）;%计算DML准则
end
[Amin,K]=min（MDL）;%寻找使DML准则最小的索引
N2=K-1;%信号源个数的估计
（2）根据
（1）中信源个数画出相应的MUSIC频率估计谱线。
%计算MUSIC谱
En=U（:
N1+1:
M）;%噪声子空间的向量组成的矩阵
NF=2048;
forn=1:
NF
Aq=exp（-j*2*pi*（n-1）/NF*（0:
M-1）'）;
Pmusic（n）=1/（Aq'*En*En'*Aq）;%music谱
Pmusic（n）=10*log10（Pmusic（n））;
end
f=-0.5:
1/NF:
0.5-1/NF;
plot（f,Pmusic）;
图7MUSIC谱估计
3.20复正弦加白噪声随机过程
同题3.18中所给。
试使用Root-MUSIC方法进行信号频率估计的仿真实验。
（要求：
信号样本数取1000，估计的自相关矩阵为8阶。
）
（1）计算正弦信号的频率估计值。
N=1000;
noise=（randn（1,N）+j*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;
fork=1:
N-M
xs（:
k）=un（k+M-1:
-1:
k）.';
end
R=xs*xs'/（N-M）;
[U,E]=svd（R）;
ev=diag（E）;
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）;
输出结果为：
最接近单位圆的两个根是：
-0.000102126882707680+1.00226420141169i
-0.000101665975924671+0.997740903255518i
上述根的归一化频率为：
0.250016217278963
0.250016217278959
（2）与MUSIC算法的估计结果比较。
En=G;
NF=2048;
forn=1:
NF
Aq=exp（-j*2*pi*（n-1）/NF*（0:
M-1）'）;
Pmusic（n）=1/（Aq'*En*En'*Aq）;
Pmusic（n）=10*log10（Pmusic（n））;
end
f=-0.5:
1/NF:
0.5-1/NF;
plot（f,Pmusic）;
图8MUSIC谱估计
3.21复正弦加白噪声随机过程
同题3.18中所给。
试使用ESPRIT算法进行信号频率估计的仿真实验，给出正弦信号频率的估计值（要求：
信号样本数取1000，估计的自相关矩阵为8阶。
）
解：
N=1000;
noise=（randn（1,N）+j*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;
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）;%计算自相关矩阵Rxx
Rxy=xs（:
1:
end-1）*xs（:
2:
end）'/（N-M-1）;%计算互相关矩阵Rxy
%—————相关矩阵的特征值分解，寻找最小特征值————————%
[U,E]=svd（Rxx）;%矩阵U是特征矢量组成的矩阵，E是对角阵，对角元素是由大到小排列的特征值
ev=diag（E）;%提取对角元素上的特征值；
emin=ev（end）;%获取最小特征值；
%——————构造矩阵对{Cxx,Cxy}————————%
Z=[zeros（M-1,1）,eye（M-1）;0,zeros（1,M-1）];%构造矩阵Z;
Cxx=Rxx-emin*eye（M）;%计算Cxx;
Cxy=Rxy-emin*Z;%计算Cxy;
%——————矩阵对{Cxx,Cxy}的广义特征值分解————————%
[U,E]=eig（Cxx,Cxy）;%广义特征值分解
z=diag（E）;%提取对角矩阵中的特征值
ph=angle（z）/（2*pi）;%求所有特征值的相位对应的归一化频率
err=abs（abs（z）-1）;%与单位元的距离
输出结果为：
单次ESPRIT算法中最接近单位圆的两个特征值为:
0.588516696017445-0.809433346373638i
-0.00139435925779186+0.999262092748743i
上述特征值的相位对应的归一化频率为：
-0.149944811600338
0.250222082900873
升级会员 