数字信号处理谱估计实验.docx
- 文档编号:7427134
- 上传时间:2023-05-11
- 格式:DOCX
- 页数:21
- 大小:27.15KB
数字信号处理谱估计实验.docx
《数字信号处理谱估计实验.docx》由会员分享,可在线阅读,更多相关《数字信号处理谱估计实验.docx(21页珍藏版)》请在冰点文库上搜索。
数字信号处理谱估计实验
实验三随机信号的功率谱估计方法
报告人:
####报告时间:
2013.1.2
一、实验目的
1.利用自相关函数法和周期图法实现对随机信号的功率谱估计。
2.观察数据长度、自相关序列长度、信噪比、窗函数、平均次数等对谱估计的分辨率、稳定性、主瓣宽度和旁瓣效应的影响。
3.学习使用FFT提高谱估计的运算速度。
4.体会非参数化功率谱估计方法的优缺点。
二、实验原理与方法
假设信号x(n为平稳随机过程,其自相关序列定义为
*([((]mExnxnmφ∆
=+(3-1
其中E表示取数学期望,*表示共轭运算。
根据定义,x(n的功率谱密度(Pω与自相关序列(mφ存在下面关系:
((jmmPmeωωφ∞-=-∞=
∑(3-2
1((2jmmPedπωπφωωπ-=⎰(3-3
然而,实际中我们很难得到准确的自相关序列(mφ,只能通过随机信号的一段样本序列来估计信号的自相关序列,进而得到信号的功率谱估计。
目前,常用的线性谱估计方法有两种,即自相关函数法和周期图法。
1.自相关函数法
假设我们已知随机信号x(n的M长的自相关序列{(mφ},利用自相关函数法可以得到x(n的功率谱估计:
*11(((LmimxiximLmφ-Λ==+-∑1
1ˆˆ((MjmmMPmeωωφ--=-+=∑(3-4利用窗函数,上式又可表达为
ˆˆ(((RjmMmP
Wmmeωωφ
∞-=-∞=∑(3-5其中,(RM
Wm为矩形窗函数,定义为1(0R
MmMWmmM<⎧=⎨≥⎩(3-6
因此,ˆ(Pω实际上是:
真正功率谱(Pω与窗函数(RM
Wm傅立叶变换的卷积。
矩形窗函数不仅降低了谱估计的分辨率,而且使谱估计产生了旁瓣。
为了降低旁瓣影响,可以采用具有较小旁瓣的窗函数,如Hamming窗,它定义为
0.540.46cos(0HMmMmWmMmMπ⎧<+⎪=⎨≥⎪⎩(3-7
这种窗函数可以有效的抑制旁瓣,但是这种方法使主瓣宽度增大,从而降低了谱估计的分辨率,这种主瓣大小和旁瓣干扰之间的矛盾在线性谱估计方法中是无法解决的。
2.周期图方法
假设已知随机信号(xn的N个样本,利用周期图方法,信号x(n的功率谱估计为
2101ˆ((NjnnPxne
Nωω--==∑(3-8
利用上述方法得到的谱估计方差与信号(xn的功率谱平方2(Pω成正比,谱估计的方差如下:
]sinsin(1[(]([22w
NNwwPwIVarxxN+⋅≈(3-9为了减小它的方差,可以将信号序列进行分段处理,然后再求各分段结果的平均,这就是平均周期图方法,即Bartlett方法。
(1Bartlett平均周期图方法
将一个随机序列(xn(0≤n≤N分成K段,每段长度为L,各段之间互不重迭,因而N=KL,可以想到,第i段的信号序列可表示为
(((1ixnxniL=+-,0,1,...,nLiK≤<=(3-10
对于每一段的周期图又可写成
210ˆ((Ljnii
nIxneωω--==∑,1,...,iK=(3-11
于是,功率谱估计定义为
11ˆˆ((K
i
iPIKωω==∑(3-12因此,对于固定的记录长度来讲,分段数K增大可使谱估计的方差减小,但是由于L的减小,相应的功率谱主瓣增宽,谱分辨率降低,显然,方差和分辨率也是矛盾的。
除了分辨率降低以外,分段处理还会引起序列的长度有限所带来的旁瓣效应。
为减小这种影响,最有效的办法是给分段序列用适当的窗函数加权,可以得到较平滑的谱估计,当然,相应的分辨率也有所下降。
(2平滑平均周期图方法
这时一种改进的Bartlett周期图方法,它特别适用于FFT直接计算功率谱估值。
将长度为N的平稳随机信号序列x(n分成K段,每段长度为L,即L=N/K。
但这里在计算周期图之前,先用窗函数(LWn给每段序列(ixn加权,K个修正的周期
图定义为
2101(((LwjniiL
nIxnWneLUωω--==∑,1,...,iK=(3-13
其中U表示窗函数序列(LWn的能量,
20
1(L
LnUWnL==∑(3-14在这种情况下,功率谱估计可按下面表达式给出:
11ˆ((KwiiPI
Kωω==∑(3-15
本实验主要是利用自相关函数法和周期图方法对下面受噪声干扰扰的正弦信号进行谱估计:
(1((iiNS
jniixnaewnωϕ+==+∑(3-16
其中NS为正弦个数,iω,iϕ和ia分别为第i个正弦信号的数字频率、
相位和幅度,iϕ随机的分布在(0,2π之间,w(n为零均值方差等于2w
σ的复高斯白噪声。
三、实验内容和分析
1.仔细阅读有关线性谱估计的内容,π根据给出的框图编制自相关函数法谱估计的程序。
运行程序,输入,0,1,0,6.0,10,1002111======waMNσϕπω选择矩形窗。
观察谱峰位置是否正确(注意:
由于窗效应可能引起谱估计的非正定。
答:
a.程序流图如3.1所示:
(程序代码见附录一
↓
↓
图表3.1实验流程框图↓↓
↓
图3.1程序流程图
b.根据代码运行结果如图3.2:
开始
输入参数:
数据长度N,自相关函数个数M,平均次数K
信号产生:
输入正弦个数Ns,每个正弦信号的数字
频率、相位和幅度,白噪声信号的方差,按照公式
(3-16产生复正弦加白噪声信号的N个采样自相关函数法
由N个x(n估计出自相
关序列(M长,并对此
自相关序列加矩形窗或
Hamming窗,利用公式(3-5计算[0,2*pi]之
间的128个功率谱抽样周期图方法输入FFT点数NF,按照公式(3-11、(3-12、(3-13计算NF点的功率谱。
结束
图3.2矩形窗进行功率谱估计的结果
分析:
由图可知,功率谱的峰值在0.6π处,与理论计算值相同。
2.观察并记录参数变化对谱估计性能的影响。
(1改变M=5,其它输入同步骤1,观察功率谱估计的主瓣宽度和旁瓣大小随自相关序列长度的变化情况。
(2选择窗函数为Hamming窗,其它输入同步骤1,观察不同的窗函数对谱估计性能的影响。
(3改变4/1πϕ=,其它输入同步骤1,观察初始相位的变化对谱估计性能的影响。
(4改变12=w
σ,其它输入同步骤1,观察信噪比变化对谱估计性能的影响。
(5改变N=10,12=w
σ,其它输入同步骤1,结合(4的内容,观察数据长度及信噪比对谱估计性能的影响。
答:
(1.M=5时(如图3.3,其他条件不变,经观察得知主瓣宽度约在(0.38-0.82,M=10时为0.48-0.70,所以主瓣宽度变宽。
同样,旁瓣宽度也明显变宽。
图3.3M=5时的矩形窗功率谱估计图
(2.当M=10其他条件不变,用Hanming窗进行功率谱估计的结果如图3.4所示:
图3.4M=10时加hanming窗的结果图
分析:
很明显,经过汉明窗处理后,主瓣宽度变大,旁瓣被抑制,但是该功率谱的分辨率比用矩形窗的结果明显降低。
(3.4/1πϕ=,其它输入同步骤1,仿真结果如图3.5所示
图3.5初始相位为45度时的功率谱估计
分析:
谱估计的峰值位置仍处于0.6π处,并且主瓣,旁瓣宽度没有明显变化,说明相位变化对谱估计的性能没有影响。
(4.改变12=w
σ,其它输入同步骤1,仿真结果如图3.6所示
图3.6噪声方差为1时的功率谱估计
分析:
高斯白噪声的方差为1时,对谱估计的峰值没有明显影响,但是影响了旁
瓣的宽度和幅值。
所以,改变信噪比对频谱的影响较小。
(5改变N=10,12=w
σ,其它输入同步骤1,其仿真结果如图3.7所示
图3.7N=10,噪声方差=1,时的功率谱估计图
分析:
结合图3.2和图3.7可以得出,当N为100时,改变信噪比对谱估计影响很小;当N较小时(实验中N=10,改变信噪比对谱估计的影响较大。
图3.7中,频谱峰值偏离了π6.0,同时旁瓣的幅值增大了许多。
3.运行程序,输入1,0,8.0,6.0,2,100212121========aaNSNϕϕπωπω,02=wσ,选择矩形窗,调整自相关序列长度M,使得两个正弦频率分量临界分辨出来,纪录此时的M值,并绘制此时的功率谱图。
同样,在加Hamming窗的情况下,记录使两个正弦频率分量临界分开的M值,并绘制此时的功率谱图。
答:
运行程序二,输入1,0,8.0,6.0,2,100212121========aaNSNϕϕπωπω,02=wσ,加矩形窗,调整M,其仿真结果如图3.8,3.9所示;加Hamming窗时分别如图3.10,3.11所示:
图3.8M=7时的俩正弦信号功率谱估计图
图3.9M=8时的俩正弦信号功率谱估计图
分析:
由图观察得知,加矩形窗时,当M=8可将两个正弦频率分量临界分辨出来。
当M=10,30时,利用汉明窗,得到的结果如图3.10,3.11所示:
图3.10M=9时汉明窗进行功率谱估计结果图
图3.11M=10时汉明窗进行功率谱估计结果图
分析:
当M=10时,可将俩个不同频率分量的信号分离。
4.运行自相关函数法谱估计程序,输入2100,10,0,1w
NMNSσ====,选择矩形窗,观察利用自相关函数法得到的白噪声信号谱估计。
改变M=3,20,观察M的变化对白噪声谱估计的影响。
答:
当M=3,10,20时结果如图3.12,3.13,3.14所示;
图3.12M=3时,矩形窗的功率谱估计
图3.13M=10时,矩形窗的功率谱估计
图3.14M=20时,矩形窗的功率谱估计
分析:
图3.12,3.13,3.14对比发现:
M越大,谱估计的各个分量变得越清晰,也就越容易分辨了。
5.根据框图,编制周期图法谱估计程序。
自程序FFT可以直接调用Matlab中的
函数。
运行程序,输入100,128,10NNFK===。
选择窗函数为矩形窗1NS=,
2
1110.6,0,0,0wwaπϕσ====,观察谱峰位置是否正确,并与步骤1结果比较。
答:
根据绘制的周期图法估计程序(见附录三,得结果如图3.15:
图3.15加窗周期图法谱估计图
分析:
对比步骤1所得的图3.2,图3.15中的谱峰位置同样也出现在0.6pi处,但
主瓣幅值明显降低,主瓣宽度增加,另外旁瓣数量减少。
6.利用周期图方法重复步骤2、3、4的内容,这里L=N/K相当于自相关函数中的M,观察周期图法谱估计和自相关函数法谱估计在分辨率和稳定性方面的差别。
答:
利用周期图法重复2,3,4的内容,(程序见附录四,这里L=N/K
相当于自相关函数法中的M,观察周期图法谱估计和自相关函数法
谱估计在分辨率和稳定性方面的差别。
N/K=10如图3.16:
N/K=3如图3.17:
N/K=20如图3.18:
图3.16N/K=10时加窗的功率谱估计
图3.17N/K=3,时加窗的功率谱估计
图3.18N/K=20,时加窗的功率谱估计
四、思考题
1.证明式(1.5可以表示为:
-1
-=0
ˆˆˆ(=2Real[(me]-(0,MjmmP
ωωφφ∑其中Real[]∙表示取实部。
证明:
10
11
1
ˆˆˆˆ((((MMjmjmjm
mMmMmP
mememe
ωωωωφ
φ
φ-----=-+=-+===+∑
∑
∑
k=-m由得:
-1
1
k00
ˆˆˆ((-(MMjk
jmmP
kemeωωωφφ
--===+∑∑ˆ-(0φ由自相关函数的对称性得:
-11
k0
ˆˆˆ(((MMjkjmmPkemeωωωφφ--===+∑∑ˆ-(0φ
即是:
10
ˆˆ(((+Mjk
jmmP
meeω
ωωφ--==∑ˆ-(0φ
即得:
-1
-=0
ˆˆˆ(=2Real[(me]-(0MjmmP
ωωφφ∑2.已知实信号(nx的N个样本(0,(1,...,(1,xxxN-,可以定义(APω的估计如下:
11
ˆˆ((,NjmA
mNPmeωωφ--=-+=∑其中||1
*0
1ˆ(((NmimxiximNφ--==+∑
试证明:
21
1ˆ((Njn
A
nPxne
N
ωω--==∑。
证明:
1||11
*1
1
1ˆˆ((((NmNNjmjmAmNimNPmexiximeNωωωφ------=-+==-+==+=
∑∑∑
1||1
*(
1
1[(](NmNjijmiimNxieximeN
ωω-----+==-++∑∑
211
*
1
1[(]((NNjnjn
jn
nnxne
xnexne
N
N
ωωω-----====∑∑#
附录一.
N=100;%样本数从0到N-1
M=10;%是矩形窗长度变量,2M-1是矩形窗长度w1=0.6*pi;
fai1=0;
a1=1;
sigma_w2=0;%方差为0
Ns=1;
I=256;%截取样本为256个
%信号产生
xn=zeros(N,1;
wn=sqrt(sigma_w2*(randn(N,1;%噪声
forn=1:
N
xn(n=a1*exp(j*(w1*n+fai1+wn(n;%定义序列end
r_xx=xcorr(xn,'unbiased';%产生自相关序列long=2*M-1;
%矩形窗
Wrect=rectwin(long;%矩形窗
p1=r_xx(N-M+1:
N+M-1.*Wrect;
p2=fft(p1,I;
p3=abs(p2(1:
I/2;%显示128个功率谱采样点figure(1,clf
stem((0:
I/2-1/(I/2,p3,'r';gridon;%归一化
title('矩形窗功率谱估计'
%hanming窗
Whanmin=hamming(long;
p11=r_xx(N-M+1:
N+M-1.*Whanmin;
p22=fft(p11,I;
p33=abs(p22(1:
I/2;
figure(2,clf
stem((0:
I/2-1/(I/2,p33,'b';gridon;
title('加hamming窗功率谱估计'
附录二.
clc;
clearall;
n=100;
m=20;
ns=0;
w1=0.6*pi;
w2=0.8*pi;
fi1=0;
fi2=0;
a1=1;
a2=1;
sigma2=1;
l=256;
sn1=zeros(n,1;
sn2=zeros(n,1;
xn=zeros(;
sn1=a1*exp(i*(w1*(1:
n'+fi1;%定义W1的信号函数
sn2=a2*exp(i*(w2*(1:
n'+fi2;%定义W2的信号函数wn=sqrt(sigma2.*(randn(n,1;%定义噪声
xn=sn1+sn2+wn;%俩个不同频率分量的信号与噪声叠加xxn=xcorr(xn,'unbiased';%产生信号
long=2*m-1;%矩形窗的定义长度为2M-1;Rectangle=window(@rectwin,long;
p1=xxn(n-m+1:
n+m-1.*Rectangle;
p2=fft(p1,l;
p3=abs(p2(1:
l/2;
long1=2*m-1;%汉明窗的定义长度为2M-1;Hammwin=window(@hamming,long1;
p11=xxn(n-m+1:
n+m-1.*Hammwin;
p22=fft(p11,l;
p33=abs(p22(1:
l/2;
figure(1;
stem((0:
l/2-1/(l/2,p3,'r';
title('矩形窗功率谱估计';gridon;figure(2;
stem((0:
l/2-1/(l/2,p33,'.';
title('汉明窗功率谱估计';gridon;附录三.
N=100;
Nf=128;
K=10;
Ns=0;
w1=0.6*pi;
fai1=0;
a1=1;
sigma_w2=0;
I=N/K;
wn=sqrt(sigma_w2*(randn(N,1;x=zeros(N,1;
fori=1:
N
x(i=a1*exp(j*(w1*i+fai1+wn(i;end
%%%加矩形窗
rectangular=window(@rectwin,I;
u=sum(rectangular.^2/I;%计算窗函数序列的能量fori=1:
K
xrectwin=x((i-1*I+1:
i*I,:
.*rectangular;xrectwindft=fft(xrectwin,Nf;
ixrec(i,:
=xrectwindft.^2/u/I;
end
pwrec=abs(sum(ixrec/K;
%%%加hamming窗
Hamming=window(@hamming,I;
u=sum(Hamming.^2/I;%计算窗函数序列能量
fori=1:
K
xhamm=x((i-1*I+1:
i*I,:
.*Hamming;
xhammdft=fft(xhamm,Nf;
ixhamm(i,:
=xhammdft.^2/Nf/I;
end
pwhamm=abs(sum(ixhamm/K;
figure(1,clf
stem((0:
Nf/2-1/(Nf/2-1,pwrec(1:
Nf/2,'r';gridon;title('加矩形窗功率谱估计'
figure(2,clf
stem((0:
Nf/2-1/(Nf/2-1,pwhamm(1:
Nf/2,'b'
title('加hamming窗功率谱估计'附录四.N=100;Nf=128;K=10;Ns=0;w1=0.6*pi;w2=0.8*pi;fai1=0;fai2=0;a1=0;a2=0;sigma_w2=1;I=N/K;wn=sqrt(sigma_w2*(randn(N,1;x=zeros(N,1;fori=1:
Nx(i=a1*exp(j*(w1*i+fai1+a2*exp(j*(w2*i+fai2+wn(i;end%%%加矩形窗rectangular=window(@rectwin,I;u=sum(rectangular.^2/I;%计算窗函数序列的能量
fori=1:
Kxrectwin=x((i-1*I+1:
i*I,:
.*rectangular;xrectwindft=fft(xrectwin,Nf;ixrec(i,:
=xrectwindft.^2/u/I;endpwrec=abs(sum(ixrec/K;%%%加hamming窗Hamming=window(@hamming,I;u=sum(Hamming.^2/I;%计算窗函数序列能量fori=1:
Kxhamm=x((i-1*I+1:
i*I,:
.*Hamming;xhammdft=fft(xhamm,Nf;ixhamm(i,:
=xhammdft.^2/Nf/I;endpwhamm=abs(sum(ixhamm/K;figure(1,clfstem((0:
Nf/2-1/(Nf/2-1,pwrec(1:
Nf/2,'r'title('加矩形窗功率谱估计'figure(2,clfstem((0:
Nf/2-1/(Nf/2-1,pwhamm(1:
Nf/2,'b'title('加hamming窗功率谱估计'
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数字信号 处理 估计 实验