自信号处理仿真题作业Word格式.docx
《自信号处理仿真题作业Word格式.docx》由会员分享，可在线阅读，更多相关《自信号处理仿真题作业Word格式.docx（18页珍藏版）》请在冰点文库上搜索。
f2=0.１７；
f３＝0.2６;
SNＲ1=30;
＞＞ＳNR2=30；
SNR3＝27;
A1＝10^（ＳNR1/20）;
A2＝10＾（SNR2/2０）;
A3=1０^（ＳＮR3/2０）;
＞sｉgnａｌ1=A1＊ｅxp（j*2*pi*f１*（0:
N－１））;
sｉgnａl2=Ａ2*exp（j*2*ｐｉ*f2*（0:
N-1））;
signal3=A3*eｘp（j*2*ｐi*ｆ3*（0：
N-１））;
un＝siｇｎal1+siｇｎal2+sｉgｎａl3+y
基于FＦT的自相关函数快速计算方法：
Ｎ=32；
Uk＝fｆt（un,2*N）；
Sｋ=（1/Ｎ）*aｂs（Uｋ）．^2;
r0=iｆft（Sk）;
r1=[r0（Ｎ+２：
2*N）,r0（1:
N）];
＞figｕre
（1）;
＞stem（real（r1））;
＞figure
（2）;
stem（ｉｍag（r1））
输出结果为：
图１基于ＦFＴ的自相关函数快速计算
实部：
虚部:
教材中式（3.1.2）估计自相关函数
r=xcorr（ｕn,N-1,'
bｉaseｄ'
）;
fｉｇure（１）;
stem（ｒｅaｌ（r））
＞figｕｒｅ
（2）;
stem（imaｇ（r））
图 2教材式（3．1．2）估计的自相关函数
（2）令信号观测样本长度N=2５6，试用BT法和周期图法估计
的功率谱，这里设BT法中所用自相关函数的单边长度M=6４。
周期图法
N=256;
noise=（rａndn（1，N）+ｊ*rａnｄn（1,N））/sqrt
（2）；
f1=0.15;
f2=0．1７;
f3=0．26;
SNR１=30；
SNＲ2=30；
SNＲ3＝2７；
A1=10^（SＮR１/2０）；
Ａ2=1０^（ＳＮR2/２０）;
Ａ3=10^（SNR3／２0）；
sｉｇｎal１＝A１*exp（j*２*pｉ＊f1*（０：
sigｎal2=Ａ2＊ｅｘp（j*2*pi*f2*（０:
sigｎal3=A３*ｅｘp（j*2*ｐi*f３*（0:
N－1））;
xn=signａl１＋sigｎaｌ2+signal3+noｉse
ＮF＝1０24;
Ｓpr=ｆftshifｔ（（1／NF）＊abs（fft（ｘn,NF））．^2）；
Spr= Spr/ｍax（Sｐr）;
t＝[-0.５:
1／NF:
０．5-１／NF];
pｌｏｔ（t,１０*lｏg10（Sｐr））;
图
３ 周期图法所得功率谱
BT法
M=64;
r＝ｘcｏｒr（xn,Ｍ，'
bｉａsed'
）；
NF＝102４;
BT＝ｆｆtｓｈift（fｆｔ（r， NF））；
BT=ＢT/ｍax（ＢＴ）；
＞＞plｏt（（-5１１:
512）/１024，10*lｏｇ10（BT））图4 BT法所得功率谱
（3）令信号观测样本长度Ｎ＝2５6,试用Leｖｉnson-Ｄurbｉn迭代算法求解AR模型的系数并估计
的功率谱,模型的阶数取为ｐ=16。
p=1６；
r0＝xcoｒr（xn,p,'
ｂiaｓed'
r=r0（ｐ＋1:
2*p+1）;
a（1,1）＝-r
（2）/r
（1）;
sigmａ
（1）=ｒ（１）-（abs（r
（2））^2）／r
（1）；
form=2：
p
k（ｍ）＝ -（r（m+1）+sｕｍ（ａ（m－1,1:
ｍ-1）.*r（m:
－１:
2）））／sigma（ｍ-1）;
a（ｍ，ｍ）=ｋ（m）;
ｆｏri＝1：
m－1
ａ（m, ｉ）=a（ｍ－１,i） +k（ｍ） *ｃｏnj（a（m-1,m-i））;
eｎd
ｓｉgma（m）=ｓiｇmａ（m-１）*（1-ａbs（k（m））^２）;
ｅnｄ
NF=１０24;
Par =sｉgmａ（ｐ）.／ffｔsｈｉfｔ（abs（ffｔ（［1,a（ｐ,:
）], ＮF））.＾2）;
Pａr＝Par/mａx（Par）;
ｐlot（（-５１1:
５12）／1024,10*log10（Par））
图5 16阶AR模型的功率谱估计
3.１８设随机过程
为
其中,
是零均值、方差为１的白噪声，
、
是相互独立并在
上服从均匀分布的随机相位。
使用ＭＶDR方法进行信号频率估计的方针实验,画出频率估计谱线，并给出正弦信号频率的估计值。
（要求：
信号样本数取1000，估计的自相关矩阵为8阶）
解:
cleaｒall;
%产生0均值,方差为1的复高斯白噪声序列v（n）
N=1000;
FS＝１;
ｎoｉse=（raｎdn（1，N）+j*ｒanｄn（1,N））/sｑｒｔ（２）；
％产生带噪声的信号样本u（n）
sｉgnal１＝eｘp（ｊ＊0．５*pｉ＊（0:
N-１）/FS+ｊ*2*ｐi＊rand）;
ｓiｇnal2=exp（-j*0.3*pi*（0：
N-1）/FS＋ｊ*2*pｉ*raｎd）;
un＝signal1＋sｉgnal2+noiｓe;
%利用un来估计自相关函数
r＝ｚｅros（1，N）;
for m=1:
N
forn=m:
r（ｍ）=r（m）＋ｕｎ（n）＊（ｕn（n-（m－1）））'
;
ｅnd
r（m）=r（m）／N;
enｄ
M=8；
％设定横向滤波器的阶数
%利用自相关矩阵和自相关函数的关系，构建自相关矩阵
R=zeroｓ（M,Ｍ）;
ｆori＝1:
M
forj=1:
iｆi<
=ｊ
R（i，j）=r（1+j-i）;
%得到了Ｍ阶的矩阵R
else
R（i,j）＝（r（1＋ｉ-ｊ））'
ｅnd
end
end
N3=2048;
%设定画图时描点的数目。
d＝１/（N3-1）；
％求画图用的横坐标的间隔。
h=zeros（1,N3）;
fori=１:
N3
h（ｉ）=-0．5+（i-１）*ｄ；
ｙ=zｅrｏs（1,Ｎ3）;
forｊ＝1:
ｗ＝h（j）＊2*ｐi;
y（j）=１０*ｌog10（aｂｓ（ｇetPMＶＤＲ（w,M，R）））;
ｅｎｄ
pｌot（h,y）;
tiｔle（＇MVＤR谱估计'
xlａbel（'
\ｏmega/2＼pｉ＇）;
ｙｌabel（＇归一化功率谱/dB'
axｉs（[－0．５ 0.５-16 1］）；
图
6MVＤR谱估计
3.１9复正弦加白噪声随机过程
同题３.18中所给。
试用ＭUSIＣ方法进行信号频率估计的仿真实验。
要求：
信号样本数取１000,估计的自相关矩阵为8阶）
（1）采用AIC和ＭＤL准则估计信号源个数；
noｉｓe=（ｒａnｄn（1，N） ＋j*raｎdn（１, N））／sqrt
（2）；
signal１ =exp（j *0.5*pi*（0:
Ｎ－1）+j*２*pi*rand）；
ｓｉｇnal2=eｘp（－ｊ*0．3＊pi *（0:
Ｎ－1）+j*2* pi *ｒａnd）;
un=signal1＋sigｎａl2+noiｓｅ；
M=８;
ｆoｒk=1：
N－M
xs（:
，k）=un（k+M-１:
-1:
k）．'
eｎd
R=ｘs*xs'
／（N－M）；
%自相关矩阵的特征值分解
[U,Ｅ]＝ｓｖd（R）;
% U是特征矢量组成的矩阵,E是对角阵,对角元素是由大到小排列的特征值
eｖ=ｄiaｇ（E）;
%提取对角元素上的特征值;
％根据AIC准则进行信号源个数的估计
fｏr k=１:
M
dｅc＝ｐrod（ev（k:
M）.^（１/（M-k+1）））;
% 计算第一项中对数的自变量的分子
nec=mｅaｎ（ev（k:
Ｍ））;
% 计算第一项中对数的自变量的分母
lnv＝（dec/nec）^（（M-k+１）*N）;
%计算第一项中对数的自变量
AIC（ｋ）＝-２*loｇ（ｌnv）+２*（k-１）*（2*M－k+１）;
%计算ＡIC准则
ｅｎd
［Ａmｉn, K]=miｎ（AIC）;
N１=Ｋ－１;
%信号源个数估计
%根据MDL准则进行信号源个数的估计
ｆｏr k=1:
ｄec=prod（ｅv（ｋ:
M）.^（1/（M-k+1）））;
nｅｃ＝mean（eｖ（k:
M））;
ｌnｖ＝（dｅc／nec）^（（M-ｋ+1）*Ｎ）;
MＤL（k）＝-log（ｌnv）+（k-1）/2*（2＊M-k+1）*log（N）；
%计算DＭL准则
eｎd
[Amin， Ｋ]=mｉｎ（MＤL）;
％寻找使ＤＭL准则最小的索引
N2=Ｋ-1;
%信号源个数的估计
（2）根据
（1）中信源个数画出相应的MUSIＣ频率估计谱线。
%计算MUSIC谱
Eｎ＝Ｕ（：
Ｎ１＋1:
M）;
%噪声子空间的向量组成的矩阵
NF=20４8;
fｏrn=1：
NＦ
Ａq=exｐ（－j*２*pi*（n-1）/NF*（0:
Ｍ－1）＇）;
Pmusic（n）＝１/（Ａq'
*En*En'
＊Ａq）；
%ｍusic谱
Ｐmusｉc（n）＝１０*loｇ１０（Pｍｕsic（ｎ））;
eｎｄ
f=-0．5:
1/NＦ:
0.５－1/NF；
plot（ｆ,Pｍusic）;
图7ＭＵSIC谱估计
３.2０复正弦加白噪声随机过程
同题3．18中所给。
试使用Roｏt－MUＳＩC方法进行信号频率估计的仿真实验。
信号样本数取１000,估计的自相关矩阵为8阶。
）
（1）计算正弦信号的频率估计值。
N＝１00０;
noｉse=（rａｎdｎ（1,Ｎ）＋j*randn（１，N））／ｓｑｒt（２）;
ｓｉgnaｌ1= ｅxp（j*0.5＊ pi*（０:
Ｎ-1）+j*2＊pi＊rand）;
signal２ =ｅxp（－j*０．3＊ｐi* （0:
N－1） +j*2*pi* raｎd）;
un=signaｌ1+sigｎal2+ｎoｉｓｅ;
M=8;
foｒk=1:
Ｎ-M
ｘs（:
，ｋ）=un（k＋Ｍ－1：
－1:
ｅnd
Ｒ=xs*xｓ'
/（Ｎ-Ｍ）;
[Ｕ，Ｅ]=sｖd（R）;
ev=ｄｉag（E）;
G＝U（:
3:
Gr=Ｇ*G'
ｃo=zerｏs（2＊M－１，1）；
ｆｏrm＝１:
cｏ（m:
m＋Ｍ-1）=co（m:
ｍ+M-1）+Gr（Ｍ:
1,m）;
eｎd
z=rｏotｓ（ｃo）；
ｐh=aｎｇｌe（ｚ）/（2＊pｉ）;
eｒｒ=abs（ａbｓ（ｚ）-1）；
输出结果为:
最接近单位圆的两个根是：
－０.27０7680+1.1６9i
-０.５92４6７１+０．９９７7４0903２５5518i
上述根的归一化频率为：
０.2563
０．2５59
（2）与MＵＳIC算法的估计结果比较。
En=G;
NF=2048;
foｒn=1:
Aq=exp（－j*2＊pｉ＊（n-1）／NF*（0:
M-1）'
Pmｕsic（n）=１/（Ａｑ＇*En*Eｎ'
*Aｑ）;
Pmusic（n）=10*log10（Pmuｓic（ｎ））;
f=-０.5:
1／ＮF:
0.５-1/NＦ;
plｏt（f,Pｍusic）；
图8 ＭUSIC谱估计
３．21复正弦加白噪声随机过程
同题3.18中所给。
试使用ESＰRIＴ算法进行信号频率估计的仿真实验，给出正弦信号频率的估计值（要求:
信号样本数取1000，估计的自相关矩阵为8阶。
N=10０0；
noｉｓｅ=（rａndn（1,N） +ｊ*raｎdｎ（1,N））/sｑｒt
（2）;
siｇｎal1＝exp（ｊ* ０．5＊ pi* （0：
Ｎ-1）+ j*2*ｐi＊ｒand）;
signal2=exp（-j*0.3＊ｐi*（0:
Ｎ－1） +j* ２*ｐi *rand）；
un=ｓiｇｎａl1+signal２+noise;
M＝8;
ｆoｒk=1:
N－M
xｓ（：
k）=un（k+M－1:
Rｘｘ=xｓ（:
，１:
eｎd-1）*xs（：
1:
end-1）'
/（N-M-1）；
%计算自相关矩阵Rxx
Rxy=xs（：
ｅｎｄ-１）*xｓ（:
２:
ｅｎd）＇/（N-M-１）;
%计算互相关矩阵Ｒxy
％—————相关矩阵的特征值分解，寻找最小特征值————————％
[Ｕ,Ｅ]=svd（Rｘx）;
％矩阵U是特征矢量组成的矩阵,E是对角阵，对角元素是由大到小排列的特征值
ev＝diag（E）;
%提取对角元素上的特征值；
ｅmiｎ=ｅv（ｅnｄ）;
%获取最小特征值；
%——————构造矩阵对｛Cｘx，Ｃxy}————————％
Z＝[zeros（Ｍ-1,1）,eye（M-1）;
0,ｚｅros（1，M－1）];
%构造矩阵Z；
Cｘx=Rｘx－ｅｍin*eｙe（M）;
%计算Cxx;
Cxy=Rxｙ-eｍｉn*Z；
%计算Ｃxy;
%——————矩阵对{Cxx,Cxy}的广义特征值分解————————%
[U,Ｅ]=eig（Cｘｘ,Ｃxｙ）;
％广义特征值分解
z＝diag（E）;
%提取对角矩阵中的特征值
ph=aｎgｌｅ（z）/（2＊pi）;
%求所有特征值的相位对应的归一化频率
eｒｒ=aｂs（ａbs（z）-１）;
%与单位元的距离
单次ESPRIT算法中最接近单位圆的两个特征值为:
0.58８5１669601７445-０.8６３8i
－０.091８６ ＋0.9９926209274８743ｉ
上述特征值的相位对应的归一化频率为：
-0．149９44811600３38
0.2５7３
升级会员 