功率谱分析

功率谱分析

【范文精选】功率谱分析

【范文大全】功率谱分析

【专家解析】功率谱分析

【优秀范文】功率谱分析

范文一:功率谱分析

随机信号处理之功率谱估计

姓名:***

学号:************ 专业:电子信息工程

一.原理

古典谱估计之相关函数法

相关法谱估计是以相关函数为媒介来计算功率谱,又叫做间接法它的理论基础是维纳-辛钦定理,其具体步骤如下:

第一步,由获得的N点数据构成的有限长序列xn(n)来估计自相关函数,即:

1N1

Rx(m)xN(n)xN(nm)

Nn0

第二步,由自相关函数的傅里叶变换求功率谱,即

ˆ(ej)Sx

M1m(M1)

ˆ(m)ejm Rx

以上两步经历了两次截断,一次是估计Rx(m) 时仅利用了x(n)的N个观测值

xN(n),这相当于对x(n)加矩形窗截断.该窗是加在数据上的,一般称为加数据窗.另

ˆ(ej)时仅利用了从-(M-1)到(M-1)的Rx(m),这相当于对Rx(m)加矩一次是估计Sxˆ(ej)分别表形窗截断,将Rx(m)截成(2M-1)长,这称为加延时窗.式中Rx(m)和Sx

j

示对Rx(m)和Sx(e)的估值.由于M≪N,使得上式的运算量不是很大,因此在FFT

问世之前,相关法是最常用的谱估计方法

古典谱估计之周期图法

首先,数据加窗得到有限长序列(www.wenku1.com)(n)直接求傅里叶变换,得频谱(www.wenku1.com)(eiω),即

1−(www.wenku1.com)

(www.wenku1.com) (www.wenku1.com) = (www.wenku1.com)− (www.wenku1.com)=0(www.wenku1.com)((www.wenku1.com))(www.wenku1.com)

然后,取频谱幅度的平方,并除以N,以此作为对x(n)真实功率谱Sx(eiω)的估计,即

1

Sx eiω =|(www.wenku1.com) (www.wenku1.com) |2

事实上,周期图法谱估计与相关法谱估计的差异只是估计自相关函数的方法不同。

参数模型法谱估计之Yule-Walker方程矩阵估计

高斯消去法是解尤勒-沃克方程的一种算法,这种方法直接求解线性方程组,运

算量较大。在解之前,先大概了解一下尤勒-沃克方程:

如前所述,P阶AR模型的系统函数为

H(z)

G1aizi

i1m

可以看出,P阶AR模型有P+1个待定系数,由上式,可得白噪声激励得到的系统输出

x(n)aix(ni)G(n)

i1

p

可以理解为,用n时刻之前的p个值的线性组合

aix(ni)

i1

p

来预测n时刻

的值x(n)预测误差为G(n).在均方误差最小准则下,组合系数

a1,a2,,ap

的选

择应使预测误差G(n)的均方值最小.经过一系列的运算,最终可以得到AR模型的正则方程

p

aiRx(mi),m1,2,pi1

Rx(m)p

aR(i)G2,m0ixi1

也就是尤勒-沃克(Yule-Walker)方程.

由P个方程可以求出P个参数(www.wenku1.com).有了参数(www.wenku1.com)(i=1,2,3,4,…,p),就可以根据自相关函数和参数(www.wenku1.com)求系统增益G。

再由公式(www.wenku1.com) (www.wenku1.com) =(www.wenku1.com) (www.wenku1.com) ∗|(www.wenku1.com)((www.wenku1.com))|(www.wenku1.com)可以得出功率谱。

参数模型法谱估计之levinson-durbin快速递推法

Levinson-durbin递推算法是解尤勒-沃克方程的快速有效的算法,这种方法利

用自相关矩阵具有的一些好的性质,是运算量大大减少,在介绍Levinson-durbin递推算法之前,先大概了解一下尤勒-沃克方程: 如前所述,P阶AR模型的系统函数为

H(z)

G1aizi

i1m

可以看出,P阶AR模型有P+1个待定系数,由上式,可得白噪声激励得到的系统输出

x(n)aix(ni)G(n)

i1

p

可以理解为,用n时刻之前的p个值的线性组合

aix(ni)

i1

p

来预测n时刻

的值x(n)预测误差为G(n).在均方误差最小准则下,组合系数

a1,a2,,ap

的选

择应使预测误差G(n)的均方值最小.经过一系列的运算,最终可以得到AR模型的正则方程

p

aiRx(mi),m1,2,pi1

Rx(m)p

aR(i)G2,m0ixi1

也就是尤勒-沃克(yule-walker)方程.

在上面方程中令p=1,得一阶AR模型的尤勒-沃克方程为

Rx(0)a1(1)Rx(1)1

Rx(1)a1(1)Rx(0)0

可以解得

a1(1)

Rx(1)Rx(0)

进而求得

R2x(1)

1Rx(0)Rx(0)[1a21(1)]

Rx(0)

在一阶的基础上进行递推,将阶次为m时的第m个系数am(m)定义为反射系数,用km表示,于是可以将计算m阶AR模型参数的Levinson-durbin递推算法表示如下:

m1

[Rx(m)am1(i)Rx(mi)]

i1

kmm1

am(i)am1(i)kmam1(mi),i1,2,3,,m1

2

mm1(1km)

式中m1,2,

2

p,0E[x(n)]Rx(0).

参数模型法谱估计之Burg算法

Burg算法不是直接估计AR模型的参数,而是先估计反射系数Km,再利用Levinson关系式求的AR模型的参数。估计反射系数所依据的准则是使前向预测误差FPE功率(www.wenku1.com)和后向预测误差BPE功率(www.wenku1.com)的平均值ε最小,即

11(www.wenku1.com)2ε= (www.wenku1.com)+(www.wenku1.com) = [(www.wenku1.com)−1 (www.wenku1.com) ]2+ [(www.wenku1.com)−1 (www.wenku1.com)−1 ] (www.wenku1.com)=(www.wenku1.com)

(www.wenku1.com)=(www.wenku1.com)

(www.wenku1.com)−1

(www.wenku1.com)−1

式中各阶前向预测误差和后向预测误差由递推公式求出,即

(www.wenku1.com) (www.wenku1.com) (www.wenku1.com)=(www.wenku1.com)−(www.wenku1.com)−1+(www.wenku1.com)−1 (www.wenku1.com) 1

(www.wenku1.com) (www.wenku1.com)0 (www.wenku1.com) =(www.wenku1.com)0(www.wenku1.com)=(www.wenku1.com)((www.wenku1.com))(www.wenku1.com)

(www.wenku1.com)

(www.wenku1.com) (www.wenku1.com) =(www.wenku1.com)−1 (www.wenku1.com) +(www.wenku1.com)−1 (www.wenku1.com)−1 (www.wenku1.com)

得到反射系数的估计公式为:

1

2 (www.wenku1.com)−(www.wenku1.com)=(www.wenku1.com)[(www.wenku1.com)−1 (www.wenku1.com) (www.wenku1.com)−1 (www.wenku1.com)−1 ]

(www.wenku1.com)=(www.wenku1.com)1 [(www.wenku1.com)22 (www.wenku1.com)− (www.wenku1.com)]+[(www.wenku1.com)−1](www.wenku1.com)=(www.wenku1.com)

(www.wenku1.com)−1

(www.wenku1.com)−1

(www.wenku1.com)

(www.wenku1.com)

上式中预测误差的求和范围表明,Burg算法所采用的数据加窗方法是协方差法,

不含有对已知数据段之外数据的人为假设。

将m阶AR模型的反射系数和m-1阶AR模型的系数代入到Levinson关系式中,可以求得AR模型其他的p-1个参数。

Levinson关系式如下:

(i)am-1(i)kmam-1(m-i),i1,2,,m-1

m阶AR模型的第m+1个参数G,

m

a

G

2

m

其中

m

是预测误差功率,可由递推公式

m



m-1

(1-km)

2

求得。

易知为进行该式的递推,必须知道0阶AR模型误差功率

j 

=

可知该式由给定序列易于求得。完成上述过程,即最终求得了表征该随机信号的AR模型的p+1个参数 。然后根据

2j H(e)(e)Sx=

即可求得该随机信号的功率谱密度。

2

Ex(n)Rx(0)

2

二.程序,图像及分析

时域波形图

i.古典谱估计之自相关法

程序: Fs=1000; N=1024; n=0:N-1;

xn=sin(2*pi*0.2*n)+0.5*sin(2*pi*0.213*n); cxn=xcorr(xn); CXk=fft(cxn,N); Pxx=abs(CXk);

index=0:round(N/2-1); k=index*Fs/N; figure(1);

plot(k, Pxx(index+1)); 图像:

图一

古典谱估计之周期图法

程序: Fs=1000; N=1024; n=0:N-1;

xn=sin(2*pi*0.2*n)+0.5*sin(2*pi*0.213*n); x1n=fft(xn); x2n=abs(x1n); x3n=x2n.*x2n; Pxx2=abs(x3n);

index=0:round(N/2-1); k=index*Fs/N; figure(2);

plot(k, Pxx2(index+1)); 图像:

图二

由图一和图二可以看出,我们利用这两种方法可以分辨出两根谱线,分别位于f=200HZ和f=213HZ附近,根据Fs=1000HZ可知,谱线位于w/pi=0.4及w/pi=0.426处,而且周期图法优于自相关法。这两种方法都是用获得的N个数据对随机过程进行功率谱估计,它隐含着对无限长数据加了一个矩形窗截断。时域与矩形窗相乘对应于频域中与矩形窗频谱卷积,故估计谱就相当于真实谱与矩形窗频谱卷积的结果。对功率谱有影响的是矩形窗频谱的幅度谱,它存在两个问题:一是主瓣不是无限窄,二是有旁瓣。由于主瓣不是无限窄,真实的功率谱与主瓣卷积后将使功率向附近频域扩散,是信号模糊,降低了谱分辨率,主瓣越宽分辨率越低。N减小时,主瓣宽度变宽,以上两图是取N=1024时得出的,此时我们假设N=128可以看出无法分辨出两根谱线。(如下图三四)由于存在旁瓣,又会产生两个后果,一是功率谱主瓣内的能量“泄露”到旁瓣使谱估计的方差增大,二是与旁瓣卷积后得到的功率谱完全属于干扰。严重情况下,强信号与旁瓣的卷积有可能大于弱信号与主瓣的卷积,使弱信号淹没在强信号的干扰中,无法检测出来。我们将w/pi=0.426处谱线的信号强度调大为w/pi=0.4处的20倍可以看出此时分辨不出两根谱线(如图五六)。可见两种方法都不是对Sx(eiω)的一致估计。于是我们利用平滑和平均的方法改进周期图法。

图三

图四

图五

图六

平滑是用一个适当的窗函数W(eiω)与算出的功率谱进行卷积,使谱线平滑。这种方法得出的谱估计是无偏的,方差也小但分辨率下降。平均是将截取的数据段(www.wenku1.com)(n)再分成L个小段,分别计算功率谱后取平均。因为L个平均的方差比随机变量单独的方差小L倍,所以当L→∞时,L个平均方差趋于0,

可以达到一致估计

的目的。

平均的程序:

N=1024;

n=0:N-1;

xn=sin(2*pi*0.2*n)+0.5*sin(2*pi*0.213*n);

Nfft=1024;

K=4;

L=N/K;

Sx=zeros(1,Nfft/2);

for k=1:K

ks=(k-1)*L+1;

ke=ks+L-1;

X=fft(x(ks:ke),Nfft);

X=abs(X).*abs(X);

for i=1:Nfft/2

Sx(i)=X(i)+Sx(i);

end

end

for i=1:Nfft/2

Sx(i)=Sx(i)/N;

end

index=0:round(Nfft/2-1);

f=index/Nfft;

plot(f,Sx(index+1));

图像:

图七

图八

显然在平均的时候(如图七),分辨率降低,且分段越多分辨率越低,如图八为分段数是14时,此时无法分辨出两根谱线。所以才采取重叠分段法,同时改变窗函数的种类观察,如下面程序即图九

N=1024;

n=0:N-1;

xn=sin(2*pi*0.2*n)+0.5*sin(2*pi*0.213*n);

Nfft=1024;

K=4;

D=fix(N/(K+1));

L=2*D;

Sx=zeros(1,Nfft/2);

wn=hamming(L);

w=wn';

for k=1:K

ks=(k-1)*D+1;

ke=ks+L-1;

xk=x(ks:ke).*w;

X=fft(xk,Nfft);

X=abs(X).*abs(X);

for i=1:Nfft/2

Sx(i)=X(i)+Sx(i);

end

end

for i=1:Nfft/2

Sx(i)=Sx(i)/N;

end

index=0:round(Nfft/2-1);

f=index/Nfft;

plot(f,Sx(index+1));

图九

在k=14时这种方法可以分辨出两根谱线,分辨率提高,如下图十

图十

参数模型法谱估计之Yule-Walker方程矩阵估计

程序:

N=1024;

n=0:N-1;

Fs=500;

x=sin(2*pi*0.2*n)+0.5*sin(2*pi*0.213*n);

for m=1:N

s=0;

for n=1:N-(m-1)

s=s+x(n)*x(n+m-1);

end

R(m)=s/N;

end

p=100;

A=ones(1,p);

for i=1:p-1

r(1,i)=A(1,i)*R(i);

end

T=toeplitz(r);

B=ones(p,1);

for i=1:p

Y(i,1)=B(i,1)*R(i+1);

end

a=-T\Y;

G=r(1,1);

for i=1:p

G=G+a(i,1)*R(i+1);

end

X=1;

K=[1 a'];

[H,w]=freqz(X,K);

Sx=G*(abs(H).*abs(H));

figure(1);

plot(w/pi,abs(Sx));

图像:

图11

12

图13

如图11可以分辨出两根谱线,位于w/pi=0.4和w/pi=0.426处,参数模型阶数的选择对求解精度产生影响。上图为P=100时的图。阶数小的时候很难识别出相邻的峰值,如图12所示是P=10时的频谱图,是估计比较平滑,相当于对实际功率谱进行平滑滤波。阶数增大时可以分辨出相邻的信号,但有可能会出现虚假峰,导致估计方差偏大。如图13所示,加了噪声之后出现虚假峰。

参数模型法谱估计之levinson-durbin快速递推法

程序:

fs=100;

N=1024 ;

n=0:N-1;

x=sin(2*pi*0.2*n)+0.5*sin(2*pi*0.213*n);

for m=1:N

R(m)=0;

end

for m=1:N

s=0;

for n=1:N-(m-1)

s=s+x(n)*x(n+m-1);

end

R(m)=s/N;

end

for M=1:N-10

a(M,M)=0;

FPE(M)=0;

P(M)=0;

end

a(1,1)=-R(2)/R(1);

P(1)=R(1)+a(1,1)*R(2);

FPE(1)=P(1)*(N+1+1)/(N-1-1);

sum=0;

for M=2:N-10

for k=1:M-1

sum=sum+a(M-1,k)*R(M-k+1);

end

a(M,M)=-(R(M+1)+sum)/P(M-1);

for k=1:M-1

a(M,k)=a(M-1,k)+a(M,M)*a(M-1,M-k);

end

P(M)=(1-(abs(a(k,k)))^2)*P(M-1);

FPE(M)=P(M)*(N+M+1)/(N-M-1);

sum=0;

end

min=FPE(1);

for M=2:N-10

if FPE(M)

min=FPE(M);

M1=M;

end

end

K=ones(1,M1+1);

for i=1:M1

K(1,i+1)=a(M1,i);

end

X=1;

[H,w]=freqz(X,K);

PSD=P(M1).*((abs(H)).^2);

plot(w/pi,abs(PSD));

图像:

图14

图15

如图14可以分辨出两根谱线,位于w/pi=0.4和w/pi=0.426

处,由矩阵估计分析

可知模型阶数对谱估计的影响,此时我们选择阶数的准则为最小预测定误差阶准则即FPE准则。阶次p的估计是使目标函数FPE(p)最小,即

FPE p =(www.wenku1.com)−(www.wenku1.com) (www.wenku1.com)随阶次增加而减小,而(www.wenku1.com)−(www.wenku1.com)N增加而增加。选择FPE(M)小于FPE(1)时的阶数为模型阶数。图15为N=512的图可以得知其分辨率并不是太高,N=128时根本无法分辨出两根谱线。由于存在谱偏移,谱分裂的问题我们采用Burg算法。 (www.wenku1.com)+(www.wenku1.com)+(www.wenku1.com)参数模型法谱估计之Burg算法

程序:

N =128;

n=1:N;

x(n) = sin(2*pi*n*0.2) +0.5*sin(2*pi*n*0.213)+0.2*randn(size(n));

p=50;

ef=zeros(p,N);

eb=zeros(p,N);

a=zeros(p,p);

ef(1,:)=x(n);

eb(1,:)=x(n);

c(1)=x*x'/N;

k(1)=0;

for m=2:p+1

x=0;

y=0;

for n=m:N

x=x+(-2)*ef(m-1,n)*eb(m-1,n-1);

y=y+(ef(m-1,n))^2+(eb(m-1,n-1))^2;

end

k(m)=x./y;

a(m,m)=k(m);

h(1)=c(1)*(1-a(1,1)^2);

for n=m:N

ef(m,n)=ef(m-1,n)+k(m).*eb(m-1,n-1);

eb(m,n)=eb(m-1,n-1)+k(m).*ef(m-1,n);

end

end

k=k(1,2:p+1);

a=a(2:p+1,2:p+1);

for m=2:p

for i=1:m-1

a(m,i)=a(m-1,i)+k(m)*a(m-1,m-i);

h(m)=h(m-1)*(1-k(m)^2);

end

end

z=a(p,:);

s=[1,z];

[H,w]=freqz(1,s);

plot(w/pi,h(p).*abs(H).^2);

图16(p=50)

图17(p=20)

图18(p=40)

图19(p=60)

图20(p=100)

图21

从图中我们可以清晰的看到Burg算法的优越性,Burg算法求解AR模型的过程是非常稳定的,而且具有很高的分辨率。当然对于Burg算法来说,P即阶数的选择是至关重要的,我们从实验结果图中可以看出,当P介于40和60之间时,得到的频谱图是较优越的,这也符合了经验定理,对于128点的频谱图分析,

P

应介于40和60之间。而当P的阶数过小的时候,会无法分辨出离的较近的两个频谱如图17,P过大频谱图会出现过多伪峰如图20,导致分辨率严重下降。在p=50时谱线分辨率较高,增加噪声功率,如图21可见会出现过多虚假峰无法分辨,信噪比较小的信号谱线丢失,说明抗干扰能力不够强。而且Burg算法中我们找不到准则来具体设定阶数P,只能根据经验。

三,实验总结

通过实验仿真可以直观地看出以下特性: 功率谱估计中古典法离散性大, 曲线粗糙, 方差较大,分辨率不太高。经典功率谱估计的分辨率反比于有效信号的长度, 但现代谱估计的分辨率可以不受此限制。这是因为对于给定的N 点有限长序列x(n), 虽然其估计出的自相关函数也是有限长的, 但是现代谱估计的一些隐含着数据和自相关函数的外推, 使其可能的长度超过给定的长度, 不象经典谱估计那样受窗函数的影响。因而现代谱的分辨率比较高, 而且现代谱线要平滑得多。

随机信号处理之功率谱估计

姓名:***

学号:************ 专业:电子信息工程

一.原理

古典谱估计之相关函数法

相关法谱估计是以相关函数为媒介来计算功率谱,又叫做间接法它的理论基础是维纳-辛钦定理,其具体步骤如下:

第一步,由获得的N点数据构成的有限长序列xn(n)来估计自相关函数,即:

1N1

Rx(m)xN(n)xN(nm)

Nn0

第二步,由自相关函数的傅里叶变换求功率谱,即

ˆ(ej)Sx

M1m(M1)

ˆ(m)ejm Rx

以上两步经历了两次截断,一次是估计Rx(m) 时仅利用了x(n)的N个观测值

xN(n),这相当于对x(n)加矩形窗截断.该窗是加在数据上的,一般称为加数据窗.另

ˆ(ej)时仅利用了从-(M-1)到(M-1)的Rx(m),这相当于对Rx(m)加矩一次是估计Sxˆ(ej)分别表形窗截断,将Rx(m)截成(2M-1)长,这称为加延时窗.式中Rx(m)和Sx

j

示对Rx(m)和Sx(e)的估值.由于M≪N,使得上式的运算量不是很大,因此在FFT

问世之前,相关法是最常用的谱估计方法

古典谱估计之周期图法

首先,数据加窗得到有限长序列

(n)直接求傅里叶变换,得频谱(eiω),即

1−

= =0()

然后,取频谱幅度的平方,并除以N,以此作为对x(n)真实功率谱Sx(eiω)的估计,即

1

Sx eiω =|

|2

事实上,周期图法谱估计与相关法谱估计的差异只是估计自相关函数的方法不同。

参数模型法谱估计之Yule-Walker方程矩阵估计

高斯消去法是解尤勒-沃克方程的一种算法,这种方法直接求解线性方程组,运

算量较大。在解之前,先大概了解一下尤勒-沃克方程:

如前所述,P阶AR模型的系统函数为

H(z)

G1aizi

i1m

可以看出,P阶AR模型有P+1个待定系数,由上式,可得白噪声激励得到的系统输出

x(n)aix(ni)G(n)

i1

p

可以理解为,用n时刻之前的p个值的线性组合

aix(ni)

i1

p

来预测n时刻

的值x(n)预测误差为G(n).在均方误差最小准则下,组合系数

a1,a2,,ap

的选

择应使预测误差G(n)的均方值最小.经过一系列的运算,最终可以得到AR模型的正则方程

p

aiRx(mi),m1,2,pi1

Rx(m)p

aR(i)G2,m0ixi1

也就是尤勒-沃克(Yule-Walker)方程.

由P个方程可以求出P个参数

.有了参数(i=1,2,3,4,…,p),就可以根据自相关函数和参数求系统增益G。

再由公式

= ∗|()|可以得出功率谱。

参数模型法谱估计之levinson-durbin快速递推法

Levinson-durbin递推算法是解尤勒-沃克方程的快速有效的算法,这种方法利

用自相关矩阵具有的一些好的性质,是运算量大大减少,在介绍Levinson-durbin递推算法之前,先大概了解一下尤勒-沃克方程: 如前所述,P阶AR模型的系统函数为

H(z)

G1aizi

i1m

可以看出,P阶AR模型有P+1个待定系数,由上式,可得白噪声激励得到的系统输出

x(n)aix(ni)G(n)

i1

p

可以理解为,用n时刻之前的p个值的线性组合

aix(ni)

i1

p

来预测n时刻

的值x(n)预测误差为G(n).在均方误差最小准则下,组合系数

a1,a2,,ap

的选

择应使预测误差G(n)的均方值最小.经过一系列的运算,最终可以得到AR模型的正则方程

p

aiRx(mi),m1,2,pi1

Rx(m)p

aR(i)G2,m0ixi1

也就是尤勒-沃克(yule-walker)方程.

在上面方程中令p=1,得一阶AR模型的尤勒-沃克方程为

Rx(0)a1(1)Rx(1)1

Rx(1)a1(1)Rx(0)0

可以解得

a1(1)

Rx(1)Rx(0)

进而求得

R2x(1)

1Rx(0)Rx(0)[1a21(1)]

Rx(0)

在一阶的基础上进行递推,将阶次为m时的第m个系数am(m)定义为反射系数,用km表示,于是可以将计算m阶AR模型参数的Levinson-durbin递推算法表示如下:

m1

[Rx(m)am1(i)Rx(mi)]

i1

kmm1

am(i)am1(i)kmam1(mi),i1,2,3,,m1

2

mm1(1km)

式中m1,2,

2

p,0E[x(n)]Rx(0).

参数模型法谱估计之Burg算法

Burg算法不是直接估计AR模型的参数,而是先估计反射系数Km,再利用Levinson关系式求的AR模型的参数。估计反射系数所依据的准则是使前向预测误差FPE功率

和后向预测误差BPE功率的平均值ε最小,即

11

2ε= + = [−1 ]2+ [−1 −1 ] =

=

−1

−1

式中各阶前向预测误差和后向预测误差由递推公式求出,即

=−1+−1 1

0 =0=()

=−1 +−1 −1

得到反射系数的估计公式为:

1

2

=[−1 −1 −1 ]

=1 [22 ]+[−1]=

−1

−1

上式中预测误差的求和范围表明,Burg算法所采用的数据加窗方法是协方差法,

不含有对已知数据段之外数据的人为假设。

将m阶AR模型的反射系数和m-1阶AR模型的系数代入到Levinson关系式中,可以求得AR模型其他的p-1个参数。

Levinson关系式如下:

(i)am-1(i)kmam-1(m-i),i1,2,,m-1

m阶AR模型的第m+1个参数G,

m

a

G

2

m

其中

m

是预测误差功率,可由递推公式

m



m-1

(1-km)

2

求得。

易知为进行该式的递推,必须知道0阶AR模型误差功率

j 

=

可知该式由给定序列易于求得。完成上述过程,即最终求得了表征该随机信号的AR模型的p+1个参数 。然后根据

2j H(e)(e)Sx=

即可求得该随机信号的功率谱密度。

2

Ex(n)Rx(0)

2

二.程序,图像及分析

时域波形图

i.古典谱估计之自相关法

程序: Fs=1000; N=1024; n=0:N-1;

xn=sin(2*pi*0.2*n)+0.5*sin(2*pi*0.213*n); cxn=xcorr(xn); CXk=fft(cxn,N); Pxx=abs(CXk);

index=0:round(N/2-1); k=index*Fs/N; figure(1);

plot(k, Pxx(index+1)); 图像:

图一

古典谱估计之周期图法

程序: Fs=1000; N=1024; n=0:N-1;

xn=sin(2*pi*0.2*n)+0.5*sin(2*pi*0.213*n); x1n=fft(xn); x2n=abs(x1n); x3n=x2n.*x2n; Pxx2=abs(x3n);

index=0:round(N/2-1); k=index*Fs/N; figure(2);

plot(k, Pxx2(index+1)); 图像:

图二

由图一和图二可以看出,我们利用这两种方法可以分辨出两根谱线,分别位于f=200HZ和f=213HZ附近,根据Fs=1000HZ可知,谱线位于w/pi=0.4及w/pi=0.426处,而且周期图法优于自相关法。这两种方法都是用获得的N个数据对随机过程进行功率谱估计,它隐含着对无限长数据加了一个矩形窗截断。时域与矩形窗相乘对应于频域中与矩形窗频谱卷积,故估计谱就相当于真实谱与矩形窗频谱卷积的结果。对功率谱有影响的是矩形窗频谱的幅度谱,它存在两个问题:一是主瓣不是无限窄,二是有旁瓣。由于主瓣不是无限窄,真实的功率谱与主瓣卷积后将使功率向附近频域扩散,是信号模糊,降低了谱分辨率,主瓣越宽分辨率越低。N减小时,主瓣宽度变宽,以上两图是取N=1024时得出的,此时我们假设N=128可以看出无法分辨出两根谱线。(如下图三四)由于存在旁瓣,又会产生两个后果,一是功率谱主瓣内的能量“泄露”到旁瓣使谱估计的方差增大,二是与旁瓣卷积后得到的功率谱完全属于干扰。严重情况下,强信号与旁瓣的卷积有可能大于弱信号与主瓣的卷积,使弱信号淹没在强信号的干扰中,无法检测出来。我们将w/pi=0.426处谱线的信号强度调大为w/pi=0.4处的20倍可以看出此时分辨不出两根谱线(如图五六)。可见两种方法都不是对Sx(eiω)的一致估计。于是我们利用平滑和平均的方法改进周期图法。

图三

图四

图五

图六

平滑是用一个适当的窗函数W(eiω)与算出的功率谱进行卷积,使谱线平滑。这种方法得出的谱估计是无偏的,方差也小但分辨率下降。平均是将截取的数据段

(n)再分成L个小段,分别计算功率谱后取平均。因为L个平均的方差比随机变量单独的方差小L倍,所以当L→∞时,L个平均方差趋于0,

可以达到一致估计

的目的。

平均的程序:

N=1024;

n=0:N-1;

xn=sin(2*pi*0.2*n)+0.5*sin(2*pi*0.213*n);

Nfft=1024;

K=4;

L=N/K;

Sx=zeros(1,Nfft/2);

for k=1:K

ks=(k-1)*L+1;

ke=ks+L-1;

X=fft(x(ks:ke),Nfft);

X=abs(X).*abs(X);

for i=1:Nfft/2

Sx(i)=X(i)+Sx(i);

end

end

for i=1:Nfft/2

Sx(i)=Sx(i)/N;

end

index=0:round(Nfft/2-1);

f=index/Nfft;

plot(f,Sx(index+1));

图像:

图七

图八

显然在平均的时候(如图七),分辨率降低,且分段越多分辨率越低,如图八为分段数是14时,此时无法分辨出两根谱线。所以才采取重叠分段法,同时改变窗函数的种类观察,如下面程序即图九

N=1024;

n=0:N-1;

xn=sin(2*pi*0.2*n)+0.5*sin(2*pi*0.213*n);

Nfft=1024;

K=4;

D=fix(N/(K+1));

L=2*D;

Sx=zeros(1,Nfft/2);

wn=hamming(L);

w=wn';

for k=1:K

ks=(k-1)*D+1;

ke=ks+L-1;

xk=x(ks:ke).*w;

X=fft(xk,Nfft);

X=abs(X).*abs(X);

for i=1:Nfft/2

Sx(i)=X(i)+Sx(i);

end

end

for i=1:Nfft/2

Sx(i)=Sx(i)/N;

end

index=0:round(Nfft/2-1);

f=index/Nfft;

plot(f,Sx(index+1));

图九

在k=14时这种方法可以分辨出两根谱线,分辨率提高,如下图十

图十

参数模型法谱估计之Yule-Walker方程矩阵估计

程序:

N=1024;

n=0:N-1;

Fs=500;

x=sin(2*pi*0.2*n)+0.5*sin(2*pi*0.213*n);

for m=1:N

s=0;

for n=1:N-(m-1)

s=s+x(n)*x(n+m-1);

end

R(m)=s/N;

end

p=100;

A=ones(1,p);

for i=1:p-1

r(1,i)=A(1,i)*R(i);

end

T=toeplitz(r);

B=ones(p,1);

for i=1:p

Y(i,1)=B(i,1)*R(i+1);

end

a=-T\Y;

G=r(1,1);

for i=1:p

G=G+a(i,1)*R(i+1);

end

X=1;

K=[1 a'];

[H,w]=freqz(X,K);

Sx=G*(abs(H).*abs(H));

figure(1);

plot(w/pi,abs(Sx));

图像:

图11

12

图13

如图11可以分辨出两根谱线,位于w/pi=0.4和w/pi=0.426处,参数模型阶数的选择对求解精度产生影响。上图为P=100时的图。阶数小的时候很难识别出相邻的峰值,如图12所示是P=10时的频谱图,是估计比较平滑,相当于对实际功率谱进行平滑滤波。阶数增大时可以分辨出相邻的信号,但有可能会出现虚假峰,导致估计方差偏大。如图13所示,加了噪声之后出现虚假峰。

参数模型法谱估计之levinson-durbin快速递推法

程序:

fs=100;

N=1024 ;

n=0:N-1;

x=sin(2*pi*0.2*n)+0.5*sin(2*pi*0.213*n);

for m=1:N

R(m)=0;

end

for m=1:N

s=0;

for n=1:N-(m-1)

s=s+x(n)*x(n+m-1);

end

R(m)=s/N;

end

for M=1:N-10

a(M,M)=0;

FPE(M)=0;

P(M)=0;

end

a(1,1)=-R(2)/R(1);

P(1)=R(1)+a(1,1)*R(2);

FPE(1)=P(1)*(N+1+1)/(N-1-1);

sum=0;

for M=2:N-10

for k=1:M-1

sum=sum+a(M-1,k)*R(M-k+1);

end

a(M,M)=-(R(M+1)+sum)/P(M-1);

for k=1:M-1

a(M,k)=a(M-1,k)+a(M,M)*a(M-1,M-k);

end

P(M)=(1-(abs(a(k,k)))^2)*P(M-1);

FPE(M)=P(M)*(N+M+1)/(N-M-1);

sum=0;

end

min=FPE(1);

for M=2:N-10

if FPE(M)

min=FPE(M);

M1=M;

end

end

K=ones(1,M1+1);

for i=1:M1

K(1,i+1)=a(M1,i);

end

X=1;

[H,w]=freqz(X,K);

PSD=P(M1).*((abs(H)).^2);

plot(w/pi,abs(PSD));

图像:

图14

图15

如图14可以分辨出两根谱线,位于w/pi=0.4和w/pi=0.426

处,由矩阵估计分析

可知模型阶数对谱估计的影响,此时我们选择阶数的准则为最小预测定误差阶准则即FPE准则。阶次p的估计是使目标函数FPE(p)最小,即

FPE p =

随阶次增加而减小,而N增加而增加。选择FPE(M)小于FPE(1)时的阶数为模型阶数。图15为N=512的图可以得知其分辨率并不是太高,N=128时根本无法分辨出两根谱线。由于存在谱偏移,谱分裂的问题我们采用Burg算法。 ++参数模型法谱估计之Burg算法

程序:

N =128;

n=1:N;

x(n) = sin(2*pi*n*0.2) +0.5*sin(2*pi*n*0.213)+0.2*randn(size(n));

p=50;

ef=zeros(p,N);

eb=zeros(p,N);

a=zeros(p,p);

ef(1,:)=x(n);

eb(1,:)=x(n);

c(1)=x*x'/N;

k(1)=0;

for m=2:p+1

x=0;

y=0;

for n=m:N

x=x+(-2)*ef(m-1,n)*eb(m-1,n-1);

y=y+(ef(m-1,n))^2+(eb(m-1,n-1))^2;

end

k(m)=x./y;

a(m,m)=k(m);

h(1)=c(1)*(1-a(1,1)^2);

for n=m:N

ef(m,n)=ef(m-1,n)+k(m).*eb(m-1,n-1);

eb(m,n)=eb(m-1,n-1)+k(m).*ef(m-1,n);

end

end

k=k(1,2:p+1);

a=a(2:p+1,2:p+1);

for m=2:p

for i=1:m-1

a(m,i)=a(m-1,i)+k(m)*a(m-1,m-i);

h(m)=h(m-1)*(1-k(m)^2);

end

end

z=a(p,:);

s=[1,z];

[H,w]=freqz(1,s);

plot(w/pi,h(p).*abs(H).^2);

图16(p=50)

图17(p=20)

图18(p=40)

图19(p=60)

图20(p=100)

图21

从图中我们可以清晰的看到Burg算法的优越性,Burg算法求解AR模型的过程是非常稳定的,而且具有很高的分辨率。当然对于Burg算法来说,P即阶数的选择是至关重要的,我们从实验结果图中可以看出,当P介于40和60之间时,得到的频谱图是较优越的,这也符合了经验定理,对于128点的频谱图分析,

P

应介于40和60之间。而当P的阶数过小的时候,会无法分辨出离的较近的两个频谱如图17,P过大频谱图会出现过多伪峰如图20,导致分辨率严重下降。在p=50时谱线分辨率较高,增加噪声功率,如图21可见会出现过多虚假峰无法分辨,信噪比较小的信号谱线丢失,说明抗干扰能力不够强。而且Burg算法中我们找不到准则来具体设定阶数P,只能根据经验。

三,实验总结

通过实验仿真可以直观地看出以下特性: 功率谱估计中古典法离散性大, 曲线粗糙, 方差较大,分辨率不太高。经典功率谱估计的分辨率反比于有效信号的长度, 但现代谱估计的分辨率可以不受此限制。这是因为对于给定的N 点有限长序列x(n), 虽然其估计出的自相关函数也是有限长的, 但是现代谱估计的一些隐含着数据和自相关函数的外推, 使其可能的长度超过给定的长度, 不象经典谱估计那样受窗函数的影响。因而现代谱的分辨率比较高, 而且现代谱线要平滑得多。

范文二:功率谱分析

功率谱分析

1.原理:

经典算法FFT:,先用FFT求出随机离散信号N点的DFT,再计算幅频特性的模平方,然后除以N,即得出该随机信号得功率谱估计。由于这种估计方法在把R(r)离散化的同时,使其功率谱周期化,故称之为“周期图法”,也称为经典谱估计方法。

AR模型正则方程与参数计算:由式(

)=

可知,

是白噪声的功

率谱(已知),要求功率谱S,只需求出H(z)

的能量谱Yule-Walker方程:p阶AR模型系统函数为H(z)=

只需求出和G(有时令G=1)。由p+1

个自相关函数( (0)到 (p))组

成的矩阵可以解出。

Levison-Durbin快速递推法:由R(0)和p个反射系数(

到)通过公式

求出

Burg算法:由预测误差递推公式

:

以及

=-

,m=1,2,,p

可得

多信号分类法(MUSIC法):求出自相关矩阵后进行特征分解,根据特征值判断

信号源个数,通过率。

=

式4.4.14,可观察出信号频

2.过程:

已知函数x(i) = 10*sin(2*pi*0.2*i)+20*sin(2*pi*0.3*i)+w(n)

其中w(n)为零均值方差为1的AWGN,采样点n=128,取p=64阶,用MATLAB分别通过4种方法估计功率谱,观察功率谱估计出

并与已知函数比较。

3.结果:

FFT

Yule-Walker方程

Levison-Durbin快速递推

Burg算法

MUSIC算法

4.分析:

FFT:

经典法进行谱估计,是有偏估计,由于卷积的运算过程会导致功率谱真实值的尖峰附近产生泄漏,相对地平滑了尖峰值,因此造成谱估计的失真。经典谱的主要缺点是频率分辨率低。这是由于周期图法在计算中把观测到的有限长的N个数据以外的数据认为是零,这显然与事实不符。

解Yule-Walker方程与Levison-Durbin快速递推方法所得结果完全一致(因为是同一种算法):

,结果很理想。

直接解Yule-Walker方程较简单, 但谱分辨率相对较差。使用Levinson- Durbin 递推可快速的求解AR 系数。

Burg算法:有一些噪声成分,但不影响

估计,功率谱更接近冲击,

,结果很理想。

Burg 算法是建立在数据基础之上的,避免了先计算自相关函数从而提高计算速度; 是较为通用的方法, 计算不太复杂, 且分辨率优于自相关法, 但对于白噪声加正弦信号有时会出现谱线分裂现象。

MUSIC

算法:

,结果很理想。

MUSIC算法是空间谱估计发展史上具有里程碑意义的算法,它实际上己经成为空间谱估计方法和理论的重要基石。其特点是测向精度高、分辨率高、能对多个来波信号同时测向及对瞬时短信号的测向:可以解决多径信号的DOA估计问题:可以用于高密度信号环境下的无线测向。如果噪声子空间大于信号子空间,MUSIC算法有更好的性能[1]。

以上结果为控制变量后所得,原始信号均为x(i) = 10*sin(2*pi*0.2*i)+20*sin(2*pi*0.3*i)+w(n),阶次均为p=64。(MUSIC算法为64点采样)

参考文献:

[1] 王娟.基于MUSIC的空间谱估计算法研究.硕士学位论文.2007

5.程序:

FFT:

clear; %清变量

n = 1:1:128; %采样点数

y = sin(2*pi*0.2*n)+sin(2*pi*0.3*n); %两正弦信号 x = awgn(y,15); %加噪(信噪比15dB) subplot(2,1,1); plot(n,x);

title('加噪的两个正弦信号');grid on; %显示加噪时域信号 f = fft(x); %fft变换 subplot(2,1,2);

plot(n/128,10*log10((abs(f).^2))); %画功率谱 xlabel('f');ylabel('功率谱dB');

title('fft经典算法求功率谱');grid on;

Yule-Walker方程:

clear; %清变量 q = 64; %阶次

for i=1:1:128 %采样点数

y(i) = 10*sin(2*pi*0.2*i)+20*sin(2*pi*0.3*i); %两正弦信号 x(i) = awgn(y(i),15); %加噪(信噪比15dB) end

subplot(2,1,1); plot(x);

title('加噪的两个正弦信号');grid on; %显示加噪时域信号 for m = 1:1:q s = 0;

for n = 1:1:128-m

s = s+x(n)*x(m+n); end

R(m) = s/128; %求自相关函数 用式2.3.10以减小方差 end

r = 0; for n = 1:1:128

r = r+x(n)*x(n); end

r = r/128; %求R(0) for i = 1:1:q

Rx(i,i) = r; %自相关矩阵对角线=R(0) end

for i = 2:1:q

for j = 1:1:i-1

Rx(i,j) = R(i-j); end end

for i = 1:1:q-1 for j = i+1:1:q

Rx(i,j) = R(j-i); end

end %求自相关矩阵Rx Ry=reshape(R,q,1); %求Ry

a = (-1)*inv(Rx)*Ry; %由Rx*a=-Ry解出系数矩阵a z(1) = 1; for i = 1:1:q

z(i+1) = a(i); %求系统函数分母多项式z end

[H,w] = freqz(1,z); %求系统函数 令G=1 subplot(2,1,2);

plot(w/(2*pi),10*log10((abs(H).^2))); %画功率谱 xlabel('f');ylabel('功率谱dB');

title('解Yule-Walker方程求功率谱');grid on;

Levison-Durbin快速递推法:

clear; %清变量 q = 64; %阶次

for i = 1:1:128 %采样点数

y(i) = 10*sin(2*pi*0.2*i)+20*sin(2*pi*0.3*i); %两正弦信号 x(i) = awgn(y(i),15); %加噪(信噪比15dB) end

subplot(2,1,1); plot(x);

title('加噪的两个正弦信号');grid on; %显示加噪时域信号 for m = 1:1:q s = 0;

for n = 1:1:128-m

s = s+x(n)*x(m+n); end

R(m) = s/128; %求自相关函数 用式2.3.10以减小方差 end

r = 0; for n = 1:1:128

r = r+x(n)*x(n); end

r = r/128; %求R(0)

k(1) = (-1)*R(1)/r; %k(1)=a(1)=-Rx(1)/Rx(0) 式3.3.15 a(1,1) = k(1);

p(1) = r*(1-k(1)*k(1)); %求预测误差功率初始值p1 for m = 2:1:q b = 0;

for i = 1:1:m-1

b = b+a(m-1,i)*R(m-i); end

k(m) = (-1)*(R(m)+b)/p(m-1); %求反射系数km 式3.3.16a a(m,m) = k(m); %am(m)=km for i = 1:1:m-1

a(m,i) = a(m-1,i)+k(m)*a(m-1,m-i); %求am(i) 式3.3.16b end

p(m) = p(m-1)*(1-k(m)*k(m)); %求预测误差功率pm 式3.3.16c end z(1) = 1; for i = 1:1:q

z(i+1) = a(q,i);

end %系统函数分母多项式z s = r;

for i = 1:1:q

s = s+a(q,i)*R(i); %求G平方 式3.3.11 end

g = sqrt(s); %求系统函数分子多项式 [H,w] = freqz(g,z); %求系统函数G subplot(2,1,2);

plot(w/(2*pi),10*log10((abs(H).^2))); %画功率谱 xlabel('f');ylabel('功率谱dB');

title('用Durbin快速递推法求功率谱');grid on;

Burg算法:

clear; %清变量 q = 64; %阶次

for i = 1:1:128 %采样点数

y(i) = 10*sin(2*pi*0.2*i)+20*sin(2*pi*0.3*i); %两正弦信号 x(i) = awgn(y(i),15); %加噪(信噪比15dB) end

subplot(2,1,1); plot(x);

title('加噪的两个正弦信号');grid on; %显示加噪时域信号 s1 = 0;

for n = 2:1:127

s1 = s1+x(n)*x(n-1); end

s2 = x(1)*x(1); for n = 2:1:127

s2 = s2+x(n)*x(n)+x(n-1)*x(n-1); end

k(1) = (-2)*s1/s2; %式3.3.25 ef(1,1) = x(1);

eb(1,1) = k(1)*x(1); for n = 2:1:128

ef(1,n) = x(n)+k(1)*x(n-1);

eb(1,n) = x(n-1)+k(1)*x(n); %式3.3.24 end

for m = 2:1:q s3 = 0;

for n = m:1:127

s3 = s3+ef(m-1,n)*eb(m-1,n-1); end s4 = 0;

for n = m:1:127

s4 = s4+ef(m-1,n)*ef(m-1,n)+eb(m-1,n-1)*eb(m-1,n-1); end

k(m) = (-2)*s3/s4; %求k(m) 式3.3.25 ef(m,1) = ef(m-1,1);

eb(m,1) = k(m)*ef(m-1,1); for n = 2:1:128

ef(m,n) = ef(m-1,n)+k(m)*eb(m-1,n-1);

eb(m,n) = eb(m-1,n-1)+k(m)*ef(m-1,n); %式3.3.24 end end

a(1,1) = k(1); for m = 2:1:q

for i = 1:1:m-1

a(m,i) = a(m-1,i)+k(m)*a(m-1,m-i); %求am(i) 式3.3.16b end

a(m,m) = k(m); %am(m)=km end z(1) = 1; for i = 1:1:q

z(i+1) = a(q,i); end

[H,w] = freqz(1,z); %求系统函数 令G=1 subplot(2,1,2);

plot(w/(2*pi),10*log10((abs(H).^2))); %画功率谱

xlabel('f');ylabel('功率谱dB');

title('用Burg算法求功率谱');grid on;

MUSIC算法:

clear; %清变量 N = 64; %采样点数 for i=1:1:N

y(i) = 10*sin(2*pi*0.2*i)+20*sin(2*pi*0.3*i); %两正弦信号 x(i) = awgn(y(i),15); %加噪(信噪比15dB) end

for m = 1:1:N-1 s = 0;

for n = 1:1:N-m

s = s+x(n)*x(m+n); end

R(m) = s/N; %求自相关函数 用式2.3.10以减小方差 end

r = 0; for n = 1:1:N

r = r+x(n)*x(n); end

r = r/N; %求R(0) for i = 1:1:N

Rx(i,i) = r; %自相关矩阵对角线=R(0) end

for i = 2:1:N

for j = 1:1:i-1

Rx(i,j) = R(i-j); end end

for i = 1:1:N-1 for j = i+1:1:N

Rx(i,j) = R(j-i); end

end %求自相关矩阵Rx [U,D,V] = svd(Rx); %矩阵的奇异值分解 M = 4; %观察D可得 S = 0;

for J = M+1:1:N

vj = V(:,J); %取第J列给vj

S = S+vj*vj'; %求和 式4.4.14中的求和 end

b = 0;

f = sqrt(-1); %f为虚数单位(相当于j) for w = 0:0.01:1 %w取值

for d = 1:1:N

e(d) = exp(f*(d-1)*w); %求e end

b = b+1;

e = reshape(e,N,1); %转置求ei

P(b) = 1/(e'*S*e); %由式4.4.14求Px end

c = 0:0.01:1

plot(c,P);

xlabel('w');ylabel('Px(w)');

title('用MUSIC算法求功率谱');grid on;

%画出Px观察可得正弦信号频率

范文三:功率谱分析

三、功率谱分析

字体 [大] [中] [小]

周期信号的功率谱为其双边幅值频谱的平方|cn|2;非周期信号的功率谱为其幅值谱密度的平方|X(ω)|2=X(ω)X*(ω)。

随机信号属于时域无限信号,其频率、幅值和相位为随机变量。因而,采用具有统计特性的功率谱估计进行谱分析

(一)自功率谱密度及其估计

各态历经随机信号的功率谱密度Sx(ω)与自相关函数Rx(τ)为傅里叶变换偶对,即

为了方便,也可用在非负频率范围内(ω>0)定义的单边功率谱密度Gx(ω)代替双边功率谱密度Sx(ω),两者之间的关系为

自功率谱估计可分为线性估计法与非线性估计法。前者以快速变换为基础,应用较早,也称为经典谱分析法; 后者是与时序模型结合的一种新方法,又称为现代谱分析方法。

1. 周期图

各态历经随机信号的均方值ψx为信号能量的时域描述。巴什瓦定理表明,信号能量的时域计算与频域计算相等,即 2

由此定义自功率谱密度及其估计为:

式中

表12-45 典型信号的自相关、频谱、概率密度

(续)

X(ω)为测试数据x(t)的傅里叶变换,X(k)为N个数据x(n)的离散傅里叶变换,由 FFT直接求出。由于X(k)具有周期函数的性质,所以称由此获得的自功率谱估计为周期图。

x′(r)的快速傅里叶变换可作为自功率谱估计的另一计算公式 自相关估计

以上两种估计都是自功率谱

Sx(ω)的有偏估计,只是偏差大小不同。

两种估计在时域对数据或对自相关估计进行截断,相当于加窗处理,致使谱估计成为真实功率谱(或称为真功率谱)与窗谱W(ω)的卷积,即

Ŝx(ω)=Sx(ω)*W(ω)

窗谱旁瓣的泄漏效应和卷积的作用使真功率谱的尖峰数值变化,邻近点的数值变大,造成谱估计的模糊与失真

以上两种估计的方差较大; 相距2π/N的各点估计值互不相关,故数据点数N越大,这些点的估计值的随机起伏越严重。

为改善谱估计的估计质量,在增大数据点数的同时,采用平均化处理和窗处理方法减小谱估计的方差。

2. 修正周期图

平均化和加窗处理可使改善的周期图估计质量提高,而且还可保留周期图便于应用FFT 计算速度快的优点。

(1) 平均周期图

把N点的长数据序列分成k段, 每段数据点数为M=N/k; 求得各段周期图Ŝx(ω)后再用平均法求得平均周期图Ŝxav(ω)。 i

平均处理使谱估计的方差减小为

当N一定时,段数K大则各段数据点数M小,谱估计的偏差大、方差小、谱平滑但频率分辨率低;若K小、M大,则偏差小、方差大。

(2) 加窗谱估计

选择适当的窗函数ω(r),对自相关估计序列x′(r)作加窗处理, 然后求谱估计

加窗谱估计平滑、方差小,故又称为平滑周期图。常用窗函数有矩形窗、三角窗、汉宁窗、海明窗和高斯窗等。

图12-7表示了分别用矩形窗、三角窗和汉宁窗分析同一数据所得结果。图a是被分析信号的真谱; 图b和c是用两种时宽的矩形窗分析的结果。可以看出,时域窗口窄,则频域分辨率低,两相邻谱线无法区分。图d和e表明三角窗和汉宁窗频率分辨率低,但旁瓣衰减快,谱的分布区域窄而边沿清晰。显然,应按信号的性质和处理要求适当选择窗函数。

图12-7 窗函数效果比较

(3) 平滑平均周期图

平滑谱估计计算要求所选用的窗函数保证求出的谱估计非负,但有的窗函数不满足此要求。如果先将数据x(n)分段,求出分段数据的加窗谱估计,然后再将各分段谱估计作平均化处理,即可满足谱估计非负的要求,又可减小估计偏差和估计方差,使谱估计质量提高。按此法计算的谱估计称为平滑平均周期图。

N个数据x(n)分成k段,k=N/M。按照快速傅里叶变换的要求,将段内数据补零使 M=2γ(γ为正整数)。求取各段的加窗谱估计

然后用平均法求平滑平均周期图

由于Px(k)计算用的数据经过加窗修正,为使谱估计为真谱的渐近无偏估计,在计算平滑平均谱估计Ŝx(k)时必须作相应的反修正, 上式中U就是反修正因子——归一化因子 i

(二)互功率谱密度

对于零均值各态历经随机信号x(n)、y(n)的时域和频域描述,与自相关估计、自功率谱估计类似。但是,由于互相关函数并非偶函数,所以互相关函数的估计量分别给作

互功率谱密度简称互谱,与互相关函数为傅里叶变换关系。互谱估计为互相关估计的离散傅里叶变换,可通快速傅里叶变换FFT求得。

当有限长数据x(n)、y(n)的数据点数为N时,则

估计方法和估计质量与自功率谱估计类似。

在互谱分析中, 常引用相干函数γ

间以及它们与系统特性的相互关系。 y(ω)分析两个各态历经随机信号x(t)、 y(t)之

若在某频率上,γxy2(ω)=0,则表示x(t)与y(t)在此频率上不相干。如果x(t)与y(t) 是统计独立的, 则在所有频率上都有γxy(ω)=0。若在所有频率上都有γxy(ω)=1, 则x(t) 与y(t)完全相干。当γxy2(ω)

相干函数的估计为 22

其中的自功率谱和互谱的估计都应当是经过频率平均或总体平均后的估计,而且还应保证x(t)、y(t)的测试数据是同时测量、具有相同分析参数(如记录长度T等),以免分析结果出现错误。

(三) 极大熵谱估计

伯格(J.P.Burg)提出的极大熵谱分析方法 (Maximum Entropy Method,简称 MEM)是一种新的功率谱估计方法。伯格在最大熵谱估计准则的提出和具体算法上有所创新。由此演变出来的算法有多种,被统称为“现代谱分析”。

1. 极大熵准则

传统谱估计方法实际上都是把无限长序列加窗截断后,由有限长序列获得功率谱估计。不论是对原始数据加窗还是对自相关函数加窗处理,其目的都是减小谱估计的方差、提高频率分辨率。然而窗处理不可避免地产生频域的“泄漏”,使功率谱失真。尽管在窗函数和处理方法上进行了许多研究,使得以周期图为基础的谱估计方法广泛应用,但谱估计的频率分辨率并不能令人满意。此外,传统谱估计方法通过增大数据点数来获得较高的谱估计精度。这样,不仅数据处理工作量大,而且对短记录数据或缓变信号等显然是无能为力的。

在谱估计计算中,对原始数据或自相关函数加窗处理,假设窗外的数据为零,而且还对窗内部分进行某些修正。这些人为措施增加了确定性因素,使原来具有的不确定性减少,在一定程度上歪曲了观测得到的信息。

极大熵谱分析方法在谱估计计算中不作加窗处理,而是采用适当的方法把由有限长数据求得的自相关估计外推。外推的原则是使相应的数据序列在外推点上取值的可能性具有最大的不确定性,亦即不对结果人为地增加任何附加信息。在数理统计学中,“熵”表征了各种随机试验的统计特性,是随机总体的平均不确定性的量度。极大熵谱分析法把熵的概念引入谱分析。上述外推原则就是要求在使随机过程的熵达到最大的条件下,确定未知的自相关函数值,外推原则的最大熵描述就是谱估计的极大熵准则。

2. 极大熵谱估计的基本原理

随机试验a有有限个不相容的结果Ai,相应的概率为P(Ai),且满足

随机试验的熵H(a)为其不确定性的量度。 则

对于连续随机变量,则熵为

式中的对数可以取10或取e为底,在比较熵的大小时并没有影响

正态随机变量x的概率密度为

则熵为

高斯随机过程的样本函数x(n),其数据点数为N,则其熵为

式中 |R|为N×N阶自相关矩阵[R]的行列式,也可用det[R]表示。[R]为托布列兹矩阵。

方差为σ的自噪声的熵为

Hε=-E[ln(p(ε1)p(ε2)…p(εN))]=Nlnσε

由于H随N增长将要发散,故用熵率h代替,熵率定义为

由于功率谱密度与自相关函数之间为傅里叶变换关系,则自相关矩阵[R]的行列式|R|的极限与功率谱密度Sx(f)的关系是:

式中 fN=2fc,为Nyquist采样频率;

fc——信号的频限,即带宽,

由此可以得到熵率为:

当由有限长随机信号测试数据计算自相关估计并外推时,为使每步估计的熵达到最大,则应有

(k≥m+1,m——数据点数)

由此得到:

式中,功率谱密度Sx(f)应受到熵率导数为零的条件约束,1/Sx(f)可用截短的傅里叶级数表示成为2m+1项有限级数之和,于是有

从而求得极大熵谱

式中 △t——采样时间间隔;

ak——自回归模型参数,ak*与ak共轭;

Pm——外推误差的方差

MEM谱与参数ak和Pm有关。典型的算法有Burg、Levinson、Akaike和Marple 算法。

3. 极大熵谱的特点

与经典谱估计方法相比,极大熵谱估计方法克服了窗处理带来的一系列缺陷。而且由于是连续谱,谱线光滑、谱峰陡峭,频率分辨力远高于经典谱估计; 而其谱分辨力与样本函数长度无关,特别适用于短数据样本、缓变过程的谱估计。

图12-8

表示短数据分析时FFT谱与MEM谱的对比。被分析信号是1Hz正弦波上叠加10%的白噪声。图a样本长度为1s时,两种谱均显示出主峰频率为1Hz;图b表示样本长度为0.57s时,MEM谱主峰频率清晰准确,而FFT谱却无法识别。

图12-9表示缓变信号的MEM谱与FFT谱的比较。被分析的信号是0.6Hz的正弦波,由图可见,MEM谱较为精确。

图12-8 短数据的MEM谱与FFT谱比较

图12-9 缓变数据的MEM谱与FFT谱

谱分析与谱估计在生产实践和科学研究中获得了日益广泛的应用。例如,在声纳系统中,为了寻找海洋水面舰艇或潜艇,需要对噪声信号进行谱分析,以提取有用信息,判断舰艇速度、方位、大小等; 对各种机电产品主体或部件作实际运行的谱分析,可检验设计的效果并提出改进设计的依据,或者进行在线监测与故障诊断,以便及时排除潜在的故障因素,保证安全运行。

范文四:频谱分析与功率谱分析

频谱分析(也称频率分析),是对动态信号在频率域内进行分析,分析的结果是以频率为坐标的各种物理量的谱线和曲线,可得到各种幅值以频率为变量的频谱函数F(ω)。频谱分析中可求得幅值谱、相位谱、功率谱和各种谱密度等等。频谱分析过程较为复杂,它是以傅里叶级数和傅里叶积分为基础的。

功率谱

频谱和功率谱有什么区别与联系?

谱是个很不严格的东西,常常指信号的Fourier变换,

是一个时间平均(time average)概念

功率谱的概念是针对功率有限信号的(能量有限信号可用能量谱分析),所表现的是单位频带内信号功率随频率的变换情况。保留频谱的幅度信息,但是丢掉了相位信息,所以频谱不同的信号其功率谱是可能相同的。有两个重要区别:

1。功率谱是随机过程的统计平均概念,平稳随机过程的功率谱是一个确定函数;而频谱是随机过程样本的Fourier变换,对于一个随机过程而言,频谱也是一个“随机过程”。(随机的频域序列)

2。功率概念和幅度概念的差别。此外,只能对宽平稳的各态历经的二阶矩过程谈功率谱,其存在性取决于二阶局是否存在并且二阶矩的Fourier变换收敛; 而频谱的存在性仅仅取决于该随机过程的该样本的Fourier变换是否收敛。

功率谱是个什么概念?它有单位吗?

随机信号是时域无限信号,不具备可积分条件,因此不能直接进行傅氏变换。一般用具有统计特性的功率谱来作为谱分析的依据。功率谱与自相关函数是一个傅氏变换对。功率谱具有单位频率的平均功率量纲。所以标准叫法是功率谱密度。通过功率谱密度函数,可以看出随机信号的能量随着频率的分布情况。像白噪声就是平行于w轴,在w轴上方的一条直线。

功率谱密度,从名字分解来看就是说,观察对象是功率,观察域是谱域,通常指频域,密度,就是指观察对象在观察域上的分布情况。一般我们讲的功率谱密度都是针对平稳随机过程的,由于平稳随机过程的样本函数一般不是绝对可积的,因此不能直接对它进行傅立叶分析。可以有三种办法来重新定义谱密度,来克服上述困难。

一是用相关函数的傅立叶变换来定义谱密度;二是用随机过程的有限时间傅立叶变换来定义谱密度;三是用平稳随机过程的谱分解来定义谱密度。三种定义方式对应于不同的用处,首先第一种方式前提是平稳随机过程不包含周期分量并且均值为零,这样才能保证相关函数在时差趋向于无穷时衰减,所以lonelystar说的不全对,光靠相关函数解决不了许多问题,要求太严格了;对于第二种方式,虽然一个平稳随机过程在无限时间上不能进行傅立叶变换,但是对于有限区间,傅立叶变换总是存在的,可以先架构有限时间区间上的变换,在对时间区间取极限,这个定义方式就是当前快速傅立叶变换(FFT)估计谱密度的依据;第三种方式是根据维纳的广义谐和分析理论:Generalized harmonic analysis, Acta Math, 55(1930),117-258,利用傅立叶-斯蒂吉斯积分,对均方连续的零均值平稳随机过程进行重构,在依靠正交性来建立的。

另外,对于非平稳随机过程,也有三种谱密度建立方法,由于字数限制,功率谱密度的单位是G的平方/频率。就是就是函数幅值的均方根值与频率之比。是对随机振动进行分析的重要参数。

功率谱密度的国际单位是什么?

如果是加速度功率谱密度,加速度的单位是m/s^2,

那么,加速度功率谱密度的单位就是(m/s^2)^2/Hz,

而Hz的单位是1/s,经过换算得到加速度功率谱密度的单位是m^2/s^3.

同理,如果是位移功率谱密度,它的单位就是m^2*s,

如果是弯矩功率谱密度,单位就是(N*m)^2*s

位移功率谱——m^2*s

速度功率谱——m^2/s

加速度功率谱——m^2/s^3

范文五:功率谱分析(Matlab)

功率谱分析

由题目内容,设采样频率fs=1000HZ,数据长度为256,模型阶数为14,f1=200,f2=300、250。

(1)用最大熵法进行谱估计

运行程序后,观察图像f1和f2相差较小时,功率谱变化更剧烈;模型的阶数越高,图像中能够获得的信息就越多,但同时计算量也就越大;增加数据长度可以获得更多的信息,提高了谱分析的分辨率,这是因为AR模型的谱估计隐含着对数据和自相关函数的外推,其长度可能会超过给定长度,分辨率不受信源信号的限制。

(2)分别用Levinson递推法和Burg法进行功率谱分析

 Levinson递推法

运行程序后,观察图像,f1和f2相差较小时,功率谱变化更剧烈;模型的阶数越高,图像中能够获得的信息就越多,但同时计算量也就越大;增加数据长度可以获得更多的信息,提高了谱分析的分辨率,但本题中信号为正弦信号加白噪声,故图像观察不明显。  Burg法

运行程序后,观察图像,f1和f2相差较小时,功率谱变化更剧烈;模型的阶数越高,图像中能够获得的信息就越多,但同时计算量也就越大;增加数据长度可以获得更多的信息,提高了谱分析的分辨率。

(3) 改变信号的相位、频率、信噪比,上述谱分析结果有何变化 如果正弦信号的频率过大,超过fs/2,会产生频率混叠现象,输入f1=600HZ,会在400HZ处产生一个波峰;降低信噪比会导致谱分辨率下降;信号起始相位的变动可导致谱线的偏移和分裂(我的图像观察不到)。

最大熵法估计

N=1024;

Nfft=256;

Fs=1000;

n=0:N-1;

t=n/Fs;

x1=sin(2*pi*200*t);

x2=sin(2*pi*300*t); %0.3

xn=x1+awgn(x1,10)+x2+awgn(x2,10);

[Pxx1,f]=pmem(xn,14,Nfft,Fs);

subplot(4,1,1)

plot(f,10*log10(Pxx1));

xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)');

title('MEM f2/fs=0.3,Nfft=256,Oder=14');

grid

N=1024;

Nfft=256;

Fs=1000;

n=0:N-1;

t=n/Fs;

x1=sin(2*pi*200*t);

x2=sin(2*pi*250*t); %0.25

xn=x1+awgn(x1,10)+x2+awgn(x2,10);

[Pxx1,f]=pmem(xn,14,Nfft,Fs);

subplot(4,1,2)

plot(f,10*log10(Pxx1));

xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('MEM f2/fs=0.25,Nfft=256,Oder=14');

grid

N=1024;

Nfft=512; %修改数据长度512

Fs=1000;

n=0:N-1;

t=n/Fs;

x1=sin(2*pi*200*t);

x2=sin(2*pi*300*t); %0.3

xn=x1+awgn(x1,10)+x2+awgn(x2,10);

[Pxx1,f]=pmem(xn,14,Nfft,Fs);

subplot(4,1,3)

plot(f,10*log10(Pxx1));

xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('MEM f2/fs=0.3,Nfft=512,Oder=14');

grid

N=1024;

Nfft=256;

Fs=1000;

n=0:N-1;

t=n/Fs;

x1=sin(2*pi*200*t);

x2=sin(2*pi*300*t); %0.3

xn=x1+awgn(x1,10)+x2+awgn(x2,10);

[Pxx1,f]=pmem(xn,24,Nfft,Fs); %修改阶数为24

subplot(4,1,4)

plot(f,10*log10(Pxx1));

xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('MEM f2/fs=0.3,Nfft=256,Oder=24');

Grid

Burg法估计

N=1024;

Nfft=256;

Fs=1000;

n=0:N-1;

t=n/Fs;

x1=sin(2*pi*200*t);

x2=sin(2*pi*300*t); %0.3

xn=x1+awgn(x1,10)+x2+awgn(x2,10);

[Pxx1,f]=pburg(xn,14,Nfft,Fs);

subplot(4,1,1)

plot(f,10*log10(Pxx1));

xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Burg f2/fs=300,Nfft=256, Oder=14');

grid

N=1024;

Nfft=256;

Fs=1000;

n=0:N-1;

t=n/Fs;

x1=sin(2*pi*200*t);

x2=sin(2*pi*250*t); %0.25

xn=x1+awgn(x1,10)+x2+awgn(x2,10);

[Pxx1,f]=pburg(xn,14,Nfft,Fs);

subplot(4,1,2)

plot(f,10*log10(Pxx1));

xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Burg f2/fs=250,Nfft=256, Oder=14');

grid

N=1024;

Nfft=512; %修改数据长度512

Fs=1000;

n=0:N-1;

t=n/Fs;

x1=sin(2*pi*200*t);

x2=sin(2*pi*300*t); %0.3

xn=x1+awgn(x1,10)+x2+awgn(x2,10);

[Pxx1,f]=pburg(xn,14,Nfft,Fs);

subplot(4,1,3)

plot(f,10*log10(Pxx1));

xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Burg f2/fs=300,Nfft=512, Oder=14');

grid

N=1024;

Nfft=256;

Fs=1000;

n=0:N-1;

t=n/Fs;

x1=sin(2*pi*200*t);

x2=sin(2*pi*300*t); %0.3

xn=x1+awgn(x1,10)+x2+awgn(x2,10);

[Pxx1,f]=pburg(xn,24,Nfft,Fs); %修改阶数为24 subplot(4,1,4)

plot(f,10*log10(Pxx1));

xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Burg f2/fs=300,Nfft=256, Oder=24');

grid

Levinson递推法

N=1024;

Nfft=256;

Fs=1000;

n=0:N-1;

t=n/Fs;

x1=sin(2*pi*200*t);

x2=sin(2*pi*300*t); %0.3

xn=x1+awgn(x1,10)+x2+awgn(x2,10);

[Pxx1,f]=pyulear(xn,14,Nfft,Fs);

%[Pxx1,f]=Levinson(xn,14,Nfft,Fs);

subplot(4,1,1)

plot(f,10*log10(Pxx1));

xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Levinson Nfft=256,f2/fs=0.3,Oder=14');

grid

N=1024;

Nfft=256;

Fs=1000;

n=0:N-1;

t=n/Fs;

x1=sin(2*pi*200*t);

x2=sin(2*pi*250*t); %0.25

xn=x1+awgn(x1,10)+x2+awgn(x2,10);

[Pxx1,f]=pyulear(xn,14,Nfft,Fs);

subplot(4,1,2)

plot(f,10*log10(Pxx1));

xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Levinson Nfft=256,f2/fs=0.25,Oder=14');

grid

N=1024;

Nfft=512; %修改数据长度512

Fs=1000;

n=0:N-1;

t=n/Fs;

x1=sin(2*pi*200*t);

x2=sin(2*pi*300*t); %0.3

xn=x1+awgn(x1,10)+x2+awgn(x2,10);

[Pxx1,f]=pyulear(xn,14,Nfft,Fs);

subplot(4,1,3)

plot(f,10*log10(Pxx1));

xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Levinson Nfft=512,f2/fs=0.3,Oder=14');

grid

N=1024;

Nfft=256;

Fs=1000;

n=0:N-1;

t=n/Fs;

x1=sin(2*pi*200*t);

x2=sin(2*pi*300*t); %0.3

xn=x1+awgn(x1,10)+x2+awgn(x2,10);

[Pxx1,f]=pyulear(xn,24,Nfft,Fs); %修改阶数为24 subplot(4,1,4)

plot(f,10*log10(Pxx1));

xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Levinson Nfft=256,f2/fs=0.3,Oder=24');

grid

最大熵法改变信号的相位、频率、信噪比

N=1024;

Nfft=256;

Fs=1000;

n=0:N-1;

t=n/Fs;

x1=sin(2*pi*200*t);

x2=sin(2*pi*300*t); %0.3

xn=x1+awgn(x1,10)+x2+awgn(x2,10);

[Pxx1,f]=pmem(xn,14,Nfft,Fs);

subplot(4,1,1)

plot(f,10*log10(Pxx1));

xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('MEM f2/fs=0.3,Nfft=256,Oder=14');

grid

N=1024;

Nfft=256;

Fs=1000;

n=0:N-1;

t=n/Fs;

x1=sin(2*pi*200*t+pi/6); %相位加了pi/6

x2=sin(2*pi*300*t); %0.3

xn=x1+awgn(x1,10)+x2+awgn(x2,10);

[Pxx1,f]=pmem(xn,14,Nfft,Fs);

subplot(4,1,2)

plot(f,10*log10(Pxx1));

xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('MEM f2/fs=0.3,Nfft=256,Oder=14,相位加pi/6'); grid

N=1024;

Nfft=256;

Fs=1000;

n=0:N-1;

t=n/Fs;

x1=sin(2*pi*200*t);

x2=sin(2*pi*300*t); %0.3

xn=x1+awgn(x1,5)+x2+awgn(x2,5); %性噪比改为5

[Pxx1,f]=pmem(xn,14,Nfft,Fs);

subplot(4,1,3)

plot(f,10*log10(Pxx1));

xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('MEM f2/fs=0.3,Nfft=256,Oder=14,性噪比=5'); grid

N=1024;

Nfft=256;

Fs=1000;

n=0:N-1;

t=n/Fs;

x1=sin(2*pi*300*t);

x2=sin(2*pi*400*t); %0.3

xn=x1+awgn(x1,10)+x2+awgn(x2,10);

[Pxx1,f]=pmem(xn,14,Nfft,Fs);

subplot(4,1,4)

plot(f,10*log10(Pxx1));

xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('MEM f1/fs=0.3,f2/fs=0.4,Nfft=256,Oder=14'); grid

Burg改变信号的相位、频率、信噪比

N=1024;

Nfft=256;

Fs=1000;

n=0:N-1;

t=n/Fs;

x1=sin(2*pi*200*t);

x2=sin(2*pi*300*t); %0.3

xn=x1+awgn(x1,10)+x2+awgn(x2,10);

[Pxx1,f]=pburg(xn,14,Nfft,Fs);

subplot(4,1,1)

plot(f,10*log10(Pxx1));

xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Burg f2/fs=300,Nfft=256, Oder=14');

grid

N=1024;

Nfft=256;

Fs=1000;

n=0:N-1;

t=n/Fs;

x1=sin(2*pi*200*t+pi/6); %相位加了pi/6

x2=sin(2*pi*300*t); %0.3

xn=x1+awgn(x1,10)+x2+awgn(x2,10);

[Pxx1,f]=pburg(xn,14,Nfft,Fs);

subplot(4,1,2)

plot(f,10*log10(Pxx1));

xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Burg f2/fs=300,Nfft=256, Oder=14,相位加pi/6'); grid

N=1024;

Nfft=256;

Fs=1000;

n=0:N-1;

t=n/Fs;

x1=sin(2*pi*300*t);

x2=sin(2*pi*400*t); %0.3

xn=x1+awgn(x1,10)+x2+awgn(x2,10);

[Pxx1,f]=pburg(xn,14,Nfft,Fs);

subplot(4,1,3)

plot(f,10*log10(Pxx1));

xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Burg f1/fs=300,f2/fs=400,Nfft=256, Oder=14'); grid

N=1024;

Nfft=256;

Fs=1000;

n=0:N-1;

t=n/Fs;

x1=sin(2*pi*200*t);

x2=sin(2*pi*300*t); %0.3

xn=x1+awgn(x1,5)+x2+awgn(x2,5); %性噪比改为5

[Pxx1,f]=pburg(xn,14,Nfft,Fs);

subplot(4,1,4)

plot(f,10*log10(Pxx1));

xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)');

title('Burg f2/fs=300,Nfft=256, Oder=14,性噪比=5'); grid

Levinson法改变信号的相位、频率、信噪比

N=1024;

Nfft=256;

Fs=1000;

n=0:N-1;

t=n/Fs;

x1=sin(2*pi*200*t);

x2=sin(2*pi*300*t); %0.3

xn=x1+awgn(x1,10)+x2+awgn(x2,10);

[Pxx1,f]=pyulear(xn,14,Nfft,Fs);

%[Pxx1,f]=Levinson(xn,14,Nfft,Fs);

subplot(4,1,1)

plot(f,10*log10(Pxx1));

xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Levinson Nfft=256,f2/fs=0.3,Oder=14');

grid

N=1024;

Nfft=256;

Fs=1000;

n=0:N-1;

t=n/Fs;

x1=sin(2*pi*200*t+pi/6); %相位加了pi/6

x2=sin(2*pi*300*t); %0.3

xn=x1+awgn(x1,10)+x2+awgn(x2,10);

[Pxx1,f]=pyulear(xn,14,Nfft,Fs);

%[Pxx1,f]=Levinson(xn,14,Nfft,Fs);

subplot(4,1,2)

plot(f,10*log10(Pxx1));

xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Levinson Nfft=256,f2/fs=0.3,Oder=14,相位加pi/6'); grid

N=1024;

Nfft=256;

Fs=1000;

n=0:N-1;

t=n/Fs;

x1=sin(2*pi*300*t);

x2=sin(2*pi*400*t); %0.3

xn=x1+awgn(x1,10)+x2+awgn(x2,10);

[Pxx1,f]=pyulear(xn,14,Nfft,Fs);

%[Pxx1,f]=Levinson(xn,14,Nfft,Fs);

subplot(4,1,3)

plot(f,10*log10(Pxx1));

xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Levinson Nfft=256,f1/fs=0.3,f2/fs=0.4,Oder=14'); grid

N=1024;

Nfft=256;

Fs=1000;

n=0:N-1;

t=n/Fs;

x1=sin(2*pi*200*t);

x2=sin(2*pi*300*t); %0.3

xn=x1+awgn(x1,5)+x2+awgn(x2,5); %性噪比位5

[Pxx1,f]=pyulear(xn,14,Nfft,Fs);

%[Pxx1,f]=Levinson(xn,14,Nfft,Fs);

subplot(4,1,4)

plot(f,10*log10(Pxx1));

xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)'); title('Levinson Nfft=256,f2/fs=0.3,Oder=14,性噪比=5'); grid

范文六:自功率谱和互功率谱

自功率谱和互功率谱

一、自功率谱及其性质:

1、 能量有限随机信号(确定性信号) x(n)

X(e

j

)

根据帕斯维尔定理得

x

n

2

(n)

12

X(e

j

2

)

d

自相关序列rxx(m)

n

x(n)x(nm)

r

xx

(m)

R

xx

(e

j

)

2

R

xx

(e

j

)X(e

j

)X(e

j

)

X(e

j

)

2、 无限能量随机信号x(n) (1)、表示式 x(n)的傅立叶变换不存在,令xN

x(n)

(n)

0

nNnN

N

r

xx

(0)E

x

2

(n)

x(n)x(n)

lim

12N1

x(n)x(n)

nN

N

=lim

N

N1

EnN2N1

x

N

(n)

x

N

(n) 

根据帕斯维尔定理,上式又=

lim

N

11E

2N12

xN(e

j

2

)

d 

又rxx(0)

2x

为平均功率

()d

2x

=

12

S

x

可得Sx()

(2)、性质 a、Sx(x)0

lim

N

E

2N11

xN(e

j

2

)

 

S

(x)limEx

N2N1

1

X

N

(e

j

2

)

0 

b、Sx()

S

x

()

由自相管序列性质rxx(m)得

r

xx

(m)

S

x

()

R

xx

()

R

xx

()

S

x

(

二、互功率谱

r

xy

(m)

R

xy

()

S

xy

()

*

r

xy

(m)Ex(n)



xy

y

(nm)



1、S

()

S

*xy

()

S

xy

()

m

Ex(n)



*

y

*

jm

(nm)e 

*

=m

Exny(nm)e

jn



=m

*

rxye

jn



=S2、S

xy

*xy

()

()

S

*yx

()

knm

S

xy

()

m

Ex(n)



y

*

(nm)



e

jm



m

Ex(km)



y

*

(k)



e

jm

m

Ex(km)



y

*

(k)



e

jm

*

Ey(k)x(km)em

*



jm



=S3、

*yx

() ()

S

xy

S

x

()

S

y

()自功率谱和互功率谱

一、自功率谱及其性质:

1、 能量有限随机信号(确定性信号) x(n)

X(e

j

)

根据帕斯维尔定理得

x

n

2

(n)

12

X(e

j

2

)

d

自相关序列rxx(m)

n

x(n)x(nm)

r

xx

(m)

R

xx

(e

j

)

2

R

xx

(e

j

)X(e

j

)X(e

j

)

X(e

j

)

2、 无限能量随机信号x(n) (1)、表示式 x(n)的傅立叶变换不存在,令xN

x(n)

(n)

0

nNnN

N

r

xx

(0)E

x

2

(n)

x(n)x(n)

lim

12N1

x(n)x(n)

nN

N

=lim

N

N1

EnN2N1

x

N

(n)

x

N

(n) 

根据帕斯维尔定理,上式又=

lim

N

11E

2N12

xN(e

j

2

)

d 

又rxx(0)

2x

为平均功率

()d

2x

=

12

S

x

可得Sx()

(2)、性质 a、Sx(x)0

lim

N

E

2N11

xN(e

j

2

)

 

S

(x)limEx

N2N1

1

X

N

(e

j

2

)

0 

b、Sx()

S

x

()

由自相管序列性质rxx(m)得

r

xx

(m)

S

x

()

R

xx

()

R

xx

()

S

x

(

二、互功率谱

r

xy

(m)

R

xy

()

S

xy

()

*

r

xy

(m)Ex(n)



xy

y

(nm)



1、S

()

S

*xy

()

S

xy

()

m

Ex(n)



*

y

*

jm

(nm)e 

*

=m

Exny(nm)e

jn



=m

*

rxye

jn



=S2、S

xy

*xy

()

()

S

*yx

()

knm

S

xy

()

m

Ex(n)



y

*

(nm)



e

jm



m

Ex(km)



y

*

(k)



e

jm

m

Ex(km)



y

*

(k)



e

jm

*

Ey(k)x(km)em

*



jm



=S3、

*yx

() ()

S

xy

S

x

()

S

y

()

范文七:话音信号的功率谱分析

清华大学学报(自然科学版)1999年第39卷第5期

CN1122223

N.39,No.5JTsinghuaUniv(Sci&Tech),1999,Vol3034

112~114

话音信号的功率谱分析3

马道钧

北京电子科技学院计算机科学与技术系,北京100070

文 摘 话音信号的振幅分布特性服从拉普拉斯分布或伽玛分布,并与方差有关。改变其方差即改变其分布特性,从而达到话音信号振幅平滑处理的目的,进而为抑制话音信号处理过程中信号频带的拓宽作准备。通过实验与统计分析论证,采用适当方式对话音信号进行处理,可以有效改变其振幅分布特性。这一“平滑”实现是十分有益的。

关键词 话音信号;;;开方电

路;乘方电路

分类号 TN918.4;TP391

率,简称共振峰。发辅音时,气流在声道中要经过不

同部位的阻碍,,尤其是前三个。而基音频率和谐波素。一般说来,男声的基音频率较低,女声的基音频率较高。基音频率通常在80~400Hz之间。

语音信号波形的典型例子如图1所示

众所周知,人类可闻声音的频率范围为16~20000Hz,而通过电话线路传递的语音信号,其频率通常需要限制在300~3400Hz之间。然而,在电话信号的处理过程中,因这些处理常常引起信号频带的拓宽,从而使经过电话线路传输后的信号产生较大的失真,而使接收端信号还原的过程无法实现。因此,研究语音信号的功率谱密度,探索语音信号进行某种预处理的方式,从而限制信号频带的拓展,对于语音信号的处理则是十分有益的。

图1 语音信号的典型特征

1 关于语音信号

人们在讲话时,信息在大脑中变换为指挥各器官动作的神经电脉冲信号,这些信号去控制发声器官产生语音波,其中包括了人们希望表达的原始信息。

简单说来,语音由元音和辅音组成。当气流从气管中发出时,它使声带产生振动,产生一个准周期的空气脉冲。这个脉冲激励某一形状的声道,产生了元音。受肌肉控制的声道的形状和大小略有不同,于是产生了不同的声音。通常称其谐振频率为共振峰频

收稿日期:1998203210

第一作者:男,1952年生,副教授

3基金项目:国家“九五”专用基金项目(95G106)

图1为男音发“Everysaltbreezecomesfrom

时的语音波形。从中可以看出,语音波形的thesea”

各段之间具有明显的停顿,而且各段波形具有突发性,即各段能量在开始时猛增,然后缓慢减小,这正是人们发音的特点。

2 语音信号的幅度概率分布特性

由语音波经电话送话器转换成电流的话音信号是时间和幅度均连续的模拟信号x(t),其信号频率集中在300~3400Hz之间,它可经过限带抽样变为无限长度的时间离散序列为

x(n)=

n=-∞

6

x(nT)∆(t-nT).

其中,T为采样周期,n=0,±1,±2,…,取其中一段有限长度序列x(n),n=0,1,2,…,N-1,其均值、均方值(平均功率)和方差为

马道钧: 话音信号的功率谱分析113

Λx=

Dx=

N

N-1

66

x(n),

2

x(n),

大的动态范围。但在较短的时间间隔内,语音信号趋于平稳。以英语为例,典型的平稳持续时间隔为100~300ms。

2

n=0

N

N-1

n=0

∆x=

N

N-1

6

[x(n)-Λx].

3 语音信号的处理与幅度概率分布的

n=0

关系

对语音信号进行某种预处理,目的是改善其幅度概率分布特性。因此,选用何种信号处理方式至关重要。

由于电话语音信号的幅度一般在-1V~1V之间,故考虑用开方电路或乘方电路对话音信号做预处理。

),由于开方电路或乘方,而不改变极,语音信号的幅度仍处1]之间。

不难证明,在[-1,1]区间,语音信号x(n)与开

(n),乘方后信号x″(n)之间的关系,方后信号x′满足

[

2

x(n)]≥[x(n)]≥[x(n)]

对于语音信号来说,当N→∞时,它接近于零均值序列。即Λx=0,∆x=Dx。然而,更有效的表达方式是研究语音信号幅度的概率密度分布。它不仅确定了过程的数字特征,还给出了波形的其它特性。如界值范围(或称为动态范围)与峰值,而这些特性对于某些话音处理系统的设计极为有用。

据有关文献介绍,在采用高性能录音设备及电话频带的条件下,话音信号的长时平均幅度概率分布特性的一阶近似可简单由拉普拉斯(Laplacian)分布(亦称双边指数型分布)或伽玛(Gamm)表示为

-PL(x)=1

x(

x3x8∆x

(1)(2)

Pr(x)=[2

3

x)1e

其理论上的拉普拉斯分布概率密度或伽玛概率

密度分布归一化后如图2所示。其中,虚线为实际测量得到的话音振幅分布。

(3)

式中[]表示仅对原始语音信号幅度的绝对值进行运算,而不改变极性。

设原始话音信号的均值和方差分别为Λx和∆x;经开方后信号的均值和方差分别为Λ′x和∆′x;经乘方后信号的均值和方差分别为Λ″不难x和∆″x。证明,当N→∞时,

Λx=

N

N-1

66

x(n)→0,x(n)→0,

2

x(n)→0.

n=0

图2 长时话音的振幅分布

Λ′x=Λ″x=

N

N-1

6

n=0N-1

如图2所示,话音信号的幅度概率分布特性具有以下特征:

1)幅度为零及其附近值的概率较大,这说明低电平的波形和无声的时间较多;

2)很大幅度出现的概率较小;

3)在二者之间的概率密度是振幅的单调函数;4)话音信号的最大幅度可达方差∆x的8倍;5)由式1和式2还可以看出,话音信号振幅分

N

n=0

即在[-1,1]区间内,经开方或乘方处理后,当N→∞时,信号仍满足零均值序列的特性。

此时,方差为

∆x=Dx=∆′x=D′x=∆″x=D″x=

显然

(4)∆′x≥∆x≥∆″x

为了进一步观察经处理后,振幅分布特性均值和方差的变化趋势,取一组数据作一些统计分析,如

N

N-1

6

[x(n)]2,[

2

x(n)],

n=0

N

N-1

6

n=0N-1

布特性与其方差有关。

与长期统计的结果不同,短期(20ms以内)的概率密度函数用简单的高斯函数就可以描述。这种短期与长期统计特性之间的差别是由语音信号波形的非平稳性引起的。这种非平稳性的定量表现就是在较长的时间间隔内,语音信号的短时方差具有较

N6

[x(n)]4.

n=0

114清华大学学报(自然科学版)1999,39(5)

(n)和x″表1所示。其中,信号S分别为x(n),x′

(n);信号幅度A的单位为V(伏);Λx为均值;∆x为方差。

由表1看出,经开方运算处理后,方差增大,而经乘方运算处理后,方差减小。这与理论分析的结果相吻合。

区间内,则得到的结论相反。换句话说,在不同区间内,采用的信号处理方式不同,方差变化亦不同。一般来说,在[0,1]区间内,经开方电路处理后,信号幅度相对向右界靠拢,使信号幅度的分布特性的方差变小;而在[-1,1]或[0,1)区间,开方电路使信号幅度峰值相对距离拉大,故而方差增大。而在此类区间内,乘方电路的处理则使信号幅度趋于零,使信号变化趋缓。

如果将信号幅度限制在[0,1]区间内,得出的结论不尽相同。由式(3),不难看出Λ′x≥Λx≥Λ″x,即经过开方电路处理以后,其均值增大;而经乘方运算处理后,均值减小。此时,经开方运算处理后,方差为

∆′x=N

N-1

4 结 论

1)话音信号的振幅分布特性虽然满足拉普拉

N

N-1

6

[x(n)-2

Λ′x]=

n=0

6

x(n)-N-1

22

Λ′x=Λx-Λ′x.

n=0

斯分布或伽玛分布,∆x有关。

,其幅度分布特性随。

,[0,1]区,乘方电路使信号幅度在1,1]区间或[0,1]区间内相对趋于幅度为零处。其作用类似于信号幅度平滑。因此,采用开方电路或乘方电路对话音信号进行处理,可以有效地改变信号幅度分布的方差和均值,即改变话音信号的振幅分布特性,使其在有效区间内相对“集中”。

3)开方电路或乘方电路的处理不同于一般的放大、滤波和限幅等,其目的是对话音信号进行平滑处理,在话音信号的持续时间内,使信号幅度的变化相对变小,因而使其有效带宽得以抑制。

4)综上所述,作为一种随机过程,话音信号振幅分布特性满足拉普拉斯分布或伽玛分布,但对其进行乘方或开方处理后,其幅度分布特性发生变化,这就为进一步分析其幅度分布特性的变化与频谱特性、频带宽度的改变提供了前提。因而这种预处理方式对话音信号处理系统是十分有益的。

而经乘方处理后,

∆″x=

N

N6

2

[x2(n)-x]n=0

N-1

6

2[x)-″x.

n=0

而原始信号的方差N-1∆x=6[x(n)-Λx]2=

N

n=0

N

N-1

6

2

x(n)-2

Λx2=Λ″x-Λ′x.

n=0

从中还很难看出方差的变化趋势,再取一组数

据作统计分析,如表2所示。

表2数据统计结果说明,在[0,1]区间内,满足

(5)∆″x≥∆x≥∆′x

比较式(4)和式(5),说明如果信号幅度限制在[0,1]区间内,则经开方运算处理后,方差减小,而经乘方运算处理后,方差增大。若信号幅度在[-1,1]

表1 在[-1,1]区间的试验数据

信号

x(n)

A

V

Λx

-0.5-0.25

0.40.16

-0.3-0.09

0.20.04

000

000

∆Λ

0.238

000

0-0.50.70.49

-0.8-0.64

0.90.81

-0.7-0.49

0.60.36

(n)x′(n)x″

0-0.7070.837-0.8940.949-0.8370.775-0.7070.632-0.5480.4470-0.25

0-0.0050.3730-0.0080.122

表2 在[0,1]区间的试验数据

信号

x(n)

A

V

Λx

0.60.7750.36

0.50.7070.25

0.30.5480.09

0.10.3160.01

000

0.5250.6370.4

∆x

0.1250.1180.133

000

0.40.6320.16

0.70.8370.49

0.90.9490.81

1.01.00.64

1.01.01.0

0.80.8941.0

(n)x′(n)x″

(下转第121页)

柴跃廷,等: 敏捷信息系统建模121

Keywords object2oriented;businessprocess

reengineering(BPR);unifiedmodelingLanguage(UML);objectmodel;behavioralmodel;systemmodel

的变化;5)通过修改对象类的操作适应业务处理功能的变化。

4 结 论

基于系统模型并配置合适的开发工具与支撑环境是实现信息系统敏捷性的主要途径,文中综合运用了面向对象的技术与方法、企业过程重组(BPR)的基本思想及UML(UnifiedModelingLanguage)的表示方法,从信息系统的静态结构与动态行为及分析、设计、实现等多角度出发,描述了敏捷信息系统的模型及其建立方法。

(上接第114页)

参 考 文 献

1JayantNS,CoxRV,McDermottBJ,etal.Analog

scramblingforspeechbasedonsequentialpermutationintimeandfrequency.SystemTechnicalJournal,1983,62(1):25~2WAD.nblingschemewhichdoes

参 考 文 献

1ChaiYueting,platform.2

LiFangyun,

Zhou

Yongchuan.development

CIMS

MIS2Oriented1(2):208~213

integrated

I:Discretetime.IEEE

heory,1979,IT225(3):261~274HJ,PiperFC.SecureSpeechCommunication.London:LondonAcademicPress,1986456

Tsinghuascienceandtechnology,柴跃廷,李芳芸,IS

MIS的快速开发平台.,6(11):19~22

葛建华.模拟语音置乱体制研究:[博士学位论文].西安:西安电子科技大学通信与电子系,1989杨保宁.模数模语音加密研究:[硕士学位论文].西安:西安电子科技大学通信与电子系,1991

文成义,赖金福,常国芩.现代电话,北京:电子工业出版社,1996

3RATIONALSoftwareration.UMLNotationVersion1.0,19974

Scholz2Reiter

B,

Stickel

E.

Business

Process

Modeling.Berlin:Springer.1996

Agileinformationsystemmodeling

CHAIYueting,ZHANGXiaodong,LIFangyun

DepartmentofAutomation,

TsinghuaUniversity,Beijing100084,China

Abstract Agileinformationsystemistheoneof

re2constructivityandfastadaptability.Thekeyofdevelopingagileinformationsystemisbuildingsystemmodelbyformalmethodology.Basedonthedeepanalysisofthecharacteristicsofenterprisebusinessactivities,boththefundamentalconstructsandtheirinter2relationshipoftheagileinformationsystemwerepresented.Intermsofthemethodologyofobject2oriented,theprincipleofBusinessProcessReengineering(BPR)andthenotationoftheUnifiedModelingLanguage(UML),themodelstructureofthiskindofsystemwereoutlined.Thesystemmodelstructurenotonlysubstantiallydifferentiatetherelativelystableandunstableaspectsofthebusinessactivities,butalsoexplicitlyreflecttherelationshipbetweenthemandintegratetheinformationofthesystemdevelopmentaswell,whichisthefoundationofthesystemreconstructionrapidly.

Analysisofpowerspectrumfor

speechsignal

MADaojun

BeijingInstituteofElectronicScienceandTechnology,

Beijing100070,China

Abstract TheamplitudedistributionofspeechsignalkeepstotheLaplaciandistributionorGammadistribution,anditisrelatedtoitsvariance.Theamplitudedistributionchangeswhthitsvariancesoastoachievesmoothprocessingofspeechsignalamplitudeandgetreadyfortherestraintofthebandwidthofspeechsignalthatisexpandedinencryptionordecryptionofspeechsignal.Manyexperimentsandstatisticalanalysisshowthatwithappropriatemethodtoprocessspeechsignal,itiseffectivetochangeitsamplitudedistribution.Smoothprocessingisveryhelpfulforthecompletionofencryptionordecryptionofspeechsignal.Keywords speechsignal;amplitudedistribution;

variance;mean;rootcircuit;squarecircuit

范文八:功率谱估计性能分析及Matlab仿真

功率谱估计性能分析及Matlab仿真

1 引言

随机信号在时域上是无限长的,在测量样本上也是无穷多的,因此随机信号的能量是无限的,应该用功率信号来描述。然而,功率信号不满足傅里叶变换的狄里克雷绝对可积的条件,因此严格意义上随机信号的傅里叶变换是不存在的。因此,要实现随机信号的频域分析,不能简单从频谱的概念出发进行研究,而是功率谱[1]。

信号的功率谱密度描述随机信号的功率在频域随频率的分布。利用给定的

N个样本数据估计一个平稳随机信号的功率谱密度叫做谱估计。谱估计方法分为两大类:经典谱估计和现代谱估计。经典功率谱估计如周期图法、自相关法等,其主要缺陷是描述功率谱波动的数字特征方差性能较差,频率分辨率低。方差性能差的原因是无法获得按功率谱密度定义中求均值和求极限的运算[2]。分辨率低的原因是在周期图法中,假定延迟窗以外的自相关函数全为0。这是不符合实际情况的,因而产生了较差的频率分辨率。而现代谱估计的目标都是旨在改善谱估计的分辨率,如自相关法和Burg法等。

2 经典功率谱估计

经典功率谱估计是截取较长的数据链中的一段作为工作区,而工作区之外的数据假设为0,这样就相当将数据加一窗函数,根据截取的N个样本数据估计出其功率谱[1]。

2.1 周期图法( Periodogram )

Schuster首先提出周期图法。周期图法是根据各态历经的随机过程功率谱的定义进行的谱估计。

取平稳随机信号x(n)的有限个观察值x(0),x(1),...,x(n1),求出其傅里叶变换

XN(e

j

)x(n)ejn

n0

N1

然后进行谱估计

1j2

S()XN(e)

N

周期图法应用比较广泛,主要是由于它与序列的频谱有直接的对应关系,并且可以采用FFT快速算法来计算。但是,这种方法需要对无限长的平稳随机序列进行截断,相当于对其加矩形窗,使之成为有限长数据。同时,这也意味着对自相关函数加三角窗,使功率谱与窗函数卷积,从而产生频谱泄露,容易使弱信号的主瓣被强信号的旁瓣所淹没,造成频谱的模糊和失真,使得谱分辨率较低[1]。

该方法基于Matlab实现的程序: clear all; load test x; N=4096;

Fn=-0.5:1/N:0.5-1/N; px=fft(x,N);

pmax=max(px);%归一化 px=px/pmax;

px=10*log10(px+0.000001); plot(Fn,fftshift(px));grid on;

图1 周期图法 N

4096

图2 周期图法 N128

说明:

(1) 本报告仿真中所采用的用于功率谱估计的数据文件来自参考文献[3]的test.dat。该数据为128点复序列(图3),由复数噪声加上四个复正弦组成。其归一化频率分别是:f10.15,f20.16,f30.252,f40.16。

图3 复序列 x(n),n[0,127]

(2) 从仿真图可以清晰看到,f1和f

2不能完全分开,仅在波形的顶部能看出

是两个频率分量;此外,当数据长度N太大时(图1),谱曲线呈现较大的起伏;当数据长度N太小时(图2),谱的分辨率又不好。据此,周期图法不满足一致性估计条件。

2.2 自相关法( BT法)

自相关法的理论基础是维纳—辛钦定理。1958年Blackman和Tukey给出了这一方法的具体实现。

对于平稳随机信号来说,其自相关函数是确定性函数,故其功率谱也是确定的。这样可由平稳随机离散信号的有限个离散值x(0),x(1),...,x(n1)求出自相关函数

1

Rx(m)

N

Nm1

n0

x(n)x(nm)

然后在(M,M)内对Rx(m)做傅里叶变换,得到功率谱

S()

mM

M

Rx(m)ejn

该方法基于Matlab实现的程序: clear all; load test x; N=4096;

Fn=-0.5:1/N:0.5-1/N; Mlag=64;

rx=xcorr(x,Mlag,'unbiased'); px=fft(rx,N);

pmax=max(px);%归一化 px=px/pmax;

px=10*log10(px+0.000001); plot(Fn,fftshift(px)); grid on;

图4 自相关法不加窗 M64

图5 自相关法不加窗 M

32

图6 自相关法使用汉明窗( Hamming )

说明:

(1) 该方法先由序列x(n)估计出自相关函数Rx(m),然后对Rx(m)进行傅里叶变换,便得到x(n)的功率谱估计。当延迟与数据长度之比很小时,可以有良好的估计精度。

(2) 图4是用自相关法(BT法)求出的功率谱,M64没有加窗;图5也是用自相关法(BT法)求出的功率谱,M32,没有加窗;图6同样是采用自相关法求出的功率谱,M32,使用了汉明窗。显然,自相关函数的延迟M越小,谱变得越平滑。

2.3 Welch法

该方法的基本原理是在对随机序列分段时,使每一段有部分重叠,然后对每一段数据用一个合适的窗函数进行平滑处理,最后对各段谱求平均。这样可得功率谱

1L

S()MULi1

M1n0

2

i

m

x

(n)ej

其中

U

(n)

n0M1

这里(n)为窗函数。

该方法基于Matlab实现的程序: clear all; load test x; N=4096;

Fn=-0.5:1/N:0.5-1/N;

xpsd=pwelch(x,hamming(33),16,N,'whole'); mmax=max(xpsd);%归一化 xpsd=xpsd/mmax;

xpsd=10*log10(xpsd+0.000001); plot(Fn,fftshift(xpsd)); grid on;

图7 Welch法 不叠合 使用汉明窗

( Hamming )

8 Welch法 叠合16点 使用汉明窗( Hamming )

图9 Welch法 叠合16点 使用矩形窗

( Boxcar )

图10 Welch法 叠合16点 使用布莱克曼窗( Blackman )

说明:

(1) 因为Welch法允许各段数据交叠,所以数据段数L会增加,使方差得到更大的改善,但是数据的交叠减小了每一段数据的不相关性,使方差的减小不会达到理论程度。另外,采用合适的窗函数可以减少信号的频谱泄露,同时也可以增加谱峰的宽度,从而提高分辨率。

(2) 图7是利用Welch法求出的周期图,共分四段,每段32点,没有叠合,使用了汉明窗;图8也是利用Welch法求出的周期图,共分四段,每段32点,,使用了汉明窗;图9是利用Welch法求出的周期图,共分四段,每段32点,交叠数为16,且使用了矩形窗;图10是利用Welch法求出的周期图,共分四段,每段32点,交叠数为16,使用了布莱克曼窗。从图中可以看出,由矩形窗处理的谱估计的主瓣宽度最窄,分辨率最好,但是其起伏性较大,所以其方差特性最差。由汉明窗和布莱克曼窗得到的谱估计的主瓣宽度最宽,因此其分辨率相对较差,但其旁瓣较小,大大改善了由矩形窗处理的谱估计旁瓣较大所产生的谱失真。因此,选择不同的窗函数其主瓣宽度不一样,造成谱估计的分辨率也不相同。

2.4 经典功率谱估计的性能比较

由以上的Matlab仿真图形和相关结果分析,我们得到了经典谱估计算法性能的直观比较:

(1) 周期图法得到的功率谱分辨率最高,但是方差性能最差,功率谱起伏剧烈,容易出现虚假谱峰。

(2) 自相关法(BT法)由于使用了平滑窗对周期图法估计的功率谱进行了平滑,因此方差性能较好,功率谱比周期图法估计的要平滑,但其分辨率比周期图法低。

(3) Welch平均周期图法是三种经典功率谱估计方法中方差性能最好的,估计的功率谱也最为平滑,但这是以分辨率的下降及偏差的增大为代价的。

综合上述讨论,我们可以对经典谱估计的算法作大致的总结[3]:

(1) 功率谱估计,不论是直接法还是间接法都可以用FFT快速计算,且物理概念明确,因而仍是目前较常用的谱估计方法。

(2) 谱的分辨率较低,它正比于2/N,N是所使用的数据长度。 (3) 方差性能不好,不是真实功率谱的一致估计,且N增大时,功率谱起伏加剧。

(4) 周期图的平滑和平均是和窗函数的使用紧密关联的,平滑和平均主要是用来改善周期图的方差性能,但往往又减小了分辨率且增加了偏差,没有一个窗函数能使估计的功率谱在方差、偏差和分辨率各个方面都得到改善,因此使用窗函数只是改进估计质量的一个技巧问题,并不能从根本上解决问题。

3 现代功率谱估计

由前一章的讨论我们可知,经典功率谱估计方法的方差性能较差,分辨率较低。而现代谱估计技术的目标都是旨在努力改善谱估计的分辨率。参数模型法是现代谱估计的主要内容,参数模型主要分为AR模型、MA模型和ARMA模型。由于AR模型具有一系列好的性能,因此是被研究最多并获得广泛应用的一种模型。本报告中现代功率谱估计的仿真基于的是AR模型。

3.1 自相关法

假定观察到得数据为x(0),x(1),...,x(n1),而对于无法观察到得区间(即,x(n)的样本假定为0,观测数据区间之外的数据为0,在均方n0和nN1)

ˆ是Toeplitz矩阵,误差意义下使得数据的预测误差功率最小。由于自相关矩阵Rp

而且又为正定的,故可利用Levinson-Durbin递归算法高效求解,得到AR模型

参数。

该方法基于Matlab实现的程序: clear all; load test x; N=4096;

fn=-0.5:1/N:0.5-1/N; xpsd=pyulear(x,20,N); pmax=max(xpsd); xpsd=xpsd/pmax;

xpsd=10*log10(xpsd+0.000001); plot(fn,fftshift(xpsd)); grid on;

图11 自相关法 p

10

图12 自相关法 p20

图13 自相关法 p30

说明:

(1) 图11、12和13是用自相关法求出的AR谱曲线,阶次p分别等于10,20和30。可以看出,在阶次较低时(图11),分辨率和检测能力均不好。当p30时,f1和f2处的两个正弦刚刚可以分开,在f3和f

4处的两个正弦也可以检出。

因此必须通过提高阶次来达到分辨出间隔较小的频率点的效果。

(2) AR模型的自相关法等效于对前向预测的误差序列前后加窗,加窗的结果是使得自相关法的分辨率降低。数据越短,分辨率越不好。

3.2 协方差法

协方差法与自相关法的区别主要在于预测误差功率求和式的上下限取得不同。由于协方差法对于观察区间0,N1外x(n)样本并未假定为0,故预测误差功率表达式中的x(nk)总是落在观察区间0,N1中,为此预测误差功率的求

ˆ是对此的半正定矩和上下限必须在p,N1之间。但由此得到的自相关矩阵Rp

阵,且不具有Toeplitz性质,故不能采用Levinson-Durbin递归算法求解,因此得到的AR模型可能不稳定。

该方法基于Matlab实现的程序: clear all; load test x;

N=4096;tn=-0.5:1/N:0.5-1/N; xpsd=pcov(x,10,N); pmax=max(xpsd); xpsd=xpsd/pmax;

xpsd=10*log10(xpsd+0.000001); plot(tn,fftshift(xpsd)); grid on;

图14 协方差法 p10

说明:

(1) 可以看到谱图在信号源频率处:f10.15,f20.16,f30.252,f40.16谱线狭窄突出,其他处谱线起伏较为剧烈。

(2) 采用协方差法对信号进行建模,能够较好地反映出信号真正的模型。

3.3修正的协方差法

AR谱估计的协方差算法基于的是最小化前向预测误差。而修正的协方差算法基于的是最小化前向和后向预测误差。这样使得它的误差功率的计算是在相对于协方差法多一倍的数据点上进行,这在观察数据长度很短的情况下,是非常有利的,但这要求信号在正反两个方向上呈现相同的特性。此外,由此得到的自相

ˆ不具有Toeplitz性质,故其正则方程不能采用Levinson-Durbin递归算关矩阵Rp

法求解。

该方法基于Matlab实现的程序: clear all; load test x; N=4096;

fn=-0.5:1/N:0.5-1/N;

xpsd=pmcov(x,10,N);

pmax=max(xpsd); xpsd=xpsd/pmax;

xpsd=10*log10(xpsd+0.000001); plot(fn,fftshift(xpsd));grid on;

图15 修正的协方差法 p10

图16 修正的协方差法p15

说明:

(1) 修正的协方差法较协方差法而言,谱估计图大致相同,但前者在信号源频率处的谱峰更加突出、尖锐,易于辨别。

(2) 通过图15和16的对比,我们可发现,在阶次较高的情况下能得到非常满意的结果。

3.4 Burg法

Burg算法是较早提出的建立在数据基础上的AR系数求解的有效算法。它基于最小化前向后向预测误差的同时满足Levinson-Durbin递归。对比其它的AR估计方法,Burg法避免了对自相关函数的计算,改而直接估计反射系数。

对于短数据的估计,Burg法求出的AR功率谱密度估计非常逼近于真值。另外,它能确保产生一个稳定的AR模型,并且能高效计算。

Burg法由于具有上述优点,所以分辨率比自相关法高,但对于混有白噪声的正弦信号,有时可能会出现谱线分裂现象。

该方法基于Matlab实现的程序: clear all; N=4096;

fn=-0.5:1/N:0.5-1/N; xpsd=pburg(x,10,N); pmax=max(xpsd); xpsd=xpsd/pmax;

xpsd=10*log10(xpsd+0.000001); plot(fn,fftshift(xpsd)); grid on;

图17 Burg法 p10

图18 Burg法 p15

说明:

(1) 图17和18是采用Burg算法对数据文件所做的功率谱估计,在阶次

p15的情况下,得到了非常满意的效果。

(2) 通过图13和17的对比,我们可明显看出Burg法的分辨率比自相关法高,在阶次较低的情况下(图17)也能较好的分辨出间隔小的频率点。

3.5 MUSIC法

刚才所讨论的自相关法、协方差法和Burg法都是基于参数建模的功率谱估计,而基于非参数建模的功率谱估计也是现代功率谱估计的重要内容。该方法是基于自相关矩阵的特征分析或者特征值分解的功率谱估计,它将相关矩阵的特征向量空间分解为信号子空间和噪声子空间,由此衍生出EV( Eigenvector )算法与MUSIC( Multiple Signal Classification )算法的信号功率谱估计[4]。其中EV谱估计与MUSIC算法谱估计都是基于噪声子空间的功率谱估计。这类方法对线谱(正弦信号的谱)最合适,对检测混有白噪声的正弦信号很有效,特别是低信噪比的情况。

MUSIC估计由下面方程给出

PMUSIC

1



eHVkVkHe

kM1

p1

1

kM1

p1

kHe

2

此处e是复正弦信号向量。 该方法基于Matlab实现的程序: clear all; load test x; N=4096;

fn=-0.5:1/N:0.5-1/N; xpsd=pmusic(x',10,N); pmax=max(xpsd); xpsd=xpsd/pmax;

xpsd=10*log10(xpsd+0.000001); for i=1:N

xxpsd(i)=xpsd(N+1-i); end

plot(fn,fftshift(xxpsd)); grid on;

图19 MUSIC法 p10

说明:

(1) 从谱图上看,MUSIC算法的谱图较为圆滑,在信号源频率处,谱线峰值明显,但峰带频率跨度较大。

(2) 通过13和19的对比,我们可明显看出MUSIC法的分辨率要好于AR模型的自相关法,在阶次较低的情况下也能较好的分辨出间隔小的频率点。

4 总结

本报告首先讨论了经典功率谱估计中三种主要方法的定义、算法及估计性能,其次讨论了现代功率谱估计中的主要内容——参数模型法谱估计,最后简要介绍了基于相关阵特征分解的谱估计算法。对所有的谱估计算法,报告均给出了Matlab程序、仿真结果及相应的分析。

参 考 文 献

[1] 王春兴.基于MATLAB实现经典功率谱估计 [J] .曲阜师范大学学报,2011,37(2) . [2] 宋宁,关华.经典功率谱估计及其仿真[J].现代电子技术,2008,31(11). [3] 胡广书.数字信号处理[M].北京:清华大学出版社,2003.

[4] 王福杰,潘宏霞.MATLAB中几种功率谱估计函数的分析比较与选择[J].计算机科学与技术,2009,27(6).

范文九:经典功率谱分析Matlab程序

一、直接法

clear;clc;close all; %清除变量;清屏;关闭当前图形窗口

Fs=1000;

t=0:1/Fs:1;

nfft=2048; %改变nfft的值可对比不同采样值时的谱估计效果

%****************生成信号、噪声**************%

x1=cos(2*pi*40*t)+3*cos(2*pi*45*t);%信号

x2=randn(size(t)); %噪声

x3=x1+x2; %信号+噪声

[Pxx,f]=periodogram(x3,window,nfft,Fs); %直接法

plot(f,10*log10(Pxx));title('直接法 nfft=2048');;

set(gca,'xlim',[1 120]); ;ylabel('Am/dB');

xlabel('Frequency/Hz');

二、间接法

Fs=1000;% 采样频率

n=0:1/Fs:1;% 产生含有噪声的序列

x1=cos(2*pi*40*n)+3*cos(2*pi*45*n);%信号

x2=randn(size(n)); %噪声

x3=x1+x2; %信号+噪声

nfft=1024;

cxn=xcorr(x3);% 计算序列的自相关函数

CXk=fft(cxn);

Pxx=abs(CXk);

index=0:round(nfft/2-1);

f=index*Fs/nfft;

plot_Pxx=10*log10(Pxx(index+1));

figure (1)

plot(f,plot_Pxx);

title('间接法 nfft=1024');ylabel('Am/dB');

set(gca,'xlim',[1 120]);

xlabel('Frequency/Hz');

三、Bartlett法

clear;clc;close all; %清除变量;清屏;关闭当前图形窗口

Fs=1000;

t=0:1/Fs:1;

nfft=1024;

%****************生成信号、噪声**************%

x1=cos(2*pi*40*t)+3*cos(2*pi*45*t);%信号

x2=randn(size(t)); %噪声

x3=x1+x2; %信号+噪声

window=hamming(512); %海明窗

noverlap=0; %数据无重叠

p=0.9; %置信概率

[Pxx,Pxxc]=psd(x3,nfft,Fs,window,noverlap,p);

index=0:round(nfft/2-1);

k=index*Fs/nfft;

plot_Pxx=10*log10(Pxx(index+1));

plot_Pxxc=10*log10(Pxxc(index+1));

figure(1)

plot(k,plot_Pxx);title('Bartlett法海明窗');;

set(gca,'xlim',[1 120]); ;ylabel('Am/dB');

xlabel('Frequency/Hz');

四、Welch法

clear;clc;close all; %清除变量;清屏;关闭当前图形窗口

Fs=1000;

t=0:1/Fs:1;

nfft=1024;

%****************生成信号、噪声**************%

x1=cos(2*pi*40*t)+3*cos(2*pi*45*t);%信号

x2=randn(size(t)); %噪声

x3=x1+x2; %信号+噪声

window=hamming(512); %海明窗

noverlap=128;

range='onesided'; %频率间隔为[0 Fs/2],只计算一半的频

[Pxx1,f]=pwelch(x3,window,noverlap,nfft,Fs,range);

plot_Pxx1=10*log10(Pxx1);

figure(1);

plot(f,plot_Pxx1);title('Welch法海明窗');ylabel('Am/dB');

set(gca,'xlim',[1 120]); xlabel('Frequency/Hz');

对所给的实测信号进行谱估计,本文采用了周期图法和Welch法。其程序如下所示:

一、周期图法

clear;clc;close all;%清除变量;清屏;关闭当前图形窗口

Fs=2000; %采样频率

load Chanel8Xia2Data.mat;

x1=Chanel8Xia2Gray;

x2=Chanel8Xia2Sti;

Mlag=length(x1);

%****************求功率谱密度**************%

[Pxx1,f]=periodogram(x1,window,length(x1),Fs);%直接法

[Pxx2,f]=periodogram(x2,window,length(x1),Fs);

plot_Pxx1=10*log10(Pxx1);

plot_Pxx2=10*log10(Pxx2);

%****************显示功率谱密度曲线**************%

figure(1);

plot(f,plot_Pxx1,'b');

axis([0,100,-50,40]);grid on;

title('直接法 ');;

set(gca,'xlim',[1 100]);

ylabel('Am/dB');

hold on;

plot(f,plot_Pxx2,'r');

axis([0,100,-50,40]);grid on;

title('直接法 ');

xlabel('Frequency/Hz');ylabel('Am/dB');

二、Welch法

clear;clc;close all; %清除变量;清屏;关闭当前图形窗口

Fs=2000;

load Chanel8Xia2Data.mat;

x1=Chanel8Xia2Gray;

x2=Chanel8Xia2Sti;

window1=hamming(1024); %海明窗

noverlap=256; %数据无重叠

range='onesided'; %频率间隔为[0 Fs/2],只计算一半的频、 %****************求功率谱密度**************%

[Pxx1,f]=pwelch(x1,window1,noverlap,length(x1),Fs,range);

[Pxx2,f]=pwelch(x2,window1,noverlap,length(x1),Fs,range);

plot_Pxx1=10*log10(Pxx1);

plot_Pxx2=10*log10(Pxx2);

%****************显示功率谱密度曲线**************%

figure(1);

plot(f,plot_Pxx1,'b');

axis([0,100,-30,30]);grid on;

title('Welch法海明窗');ylabel('Am/dB');

set(gca,'xlim',[1 100]);

hold on;

plot(f,plot_Pxx2,'r');

title('Welch法海明窗');ylabel('Am/dB');

xlabel('Frequency/Hz');

范文十:功率谱估计性能分析及其MATLAB实现

功率谱估计性能分析及其MATLAB实现

一、 经典功率谱估计分类简介

1. 间接法

根据维纳-辛钦定理,1958年Blackman和Turkey给出了这一方法的具体实现,即先由N个观察值xN(n),估计出自相关函数rx(m),求自相关函数傅里叶变换,以此变换结果作为对功率谱Px(w)的估计。

2. 直接法

直接法功率谱估计是间接法功率谱估计的一个特例,又称为周期图法,它是把随机信号的N个观察值xN(n)直接进行傅里叶变换,得到XN(ejw),然后取其幅值的平方,再除以N,作为对功率谱Px(w)的估计。

3. 改进的周期图法

将N点的观察值分成L个数据段,每段的数据为M,然后计算L个数据段的周期图的平均 per(w)的方差Pper(w),作为功率谱的估计,以此来改善用N点观察数据直接计算的周期图P

特性。根据分段方法的不同,又可以分为Welch法和Bartlett法。

Welch法

所分的数据段可以互相重叠,选用的数据窗可以是任意窗。

Bartlett法

所分的数据段互不重叠,选用的数据窗是矩形窗。

二、 经典功率谱估计的性能比较

1. 仿真结果

为了比较经典功率谱估计的性能,本文采用的信号是高斯白噪声加两个正弦信号,采样率Fs=1000Hz,两个正弦信号的频率分别为f1=200Hz,f2=210Hz。所用数据长度N=400. 仿真结果如下:

Figure1 经典功率谱估计的仿真结果

Figure1(a)示出了待估计信号的时域波形;

Figure2(b)示出了用该数据段直接求出的周期图,所用的数据窗为矩形窗;

Figure2(c)是用BT法(间接法)求出的功率谱曲线,对自相关函数用的平滑窗为矩形窗,长度M=128,数据没有加窗;

Figure2(d)是用BT法(间接法)求出的功率谱曲线,对自相关函数用的平滑窗为Hamming窗,长度M=64,数据没有加窗;

Figure2(e)是用Welch平均法求出的功率谱曲线,每段数据的长度为64点,重叠32点,使用的Hamming窗;

Figure2(f)是用Welch平均法求出的功率谱曲线,每段数据的长度为100点,重叠48点,使用的Hamming窗;

2. 性能比较

1) 直接法得到的功率谱分辨率最高,但是方差性能最差,功率谱起伏剧烈,容易出现

虚假谱峰;

2) 间接法由于使用了平滑窗对直接法估计的功率谱进行了平滑,因此方差性能比直接

法好,功率谱比直接法估计的要平滑,但其分辨率比直接法低。

3) Welch平均周期图法是三种经典功率谱估计方法中方差性能最好的,估计的功率谱

也最为平滑,但这是以分辨率的下降及偏差的增大为代价的。

3. 关于经典功率谱估计的总结

1) 功率谱估计,不论是直接法还是间接法都可以用FFT快速计算,且物理概念明确,因而

仍是目前较常用的谱估计方法。

2) 谱的分辨率较低,它正比于2π/N,N是所使用的数据长度。

3) 方差性能不好,不是真实功率谱的一致估计,且N增大时,功率谱起伏加剧。

4) 周期图的平滑和平均是和窗函数的使用紧密关联的,平滑和平均主要是用来改善周期图

的方差性能,但往往又减小了分辨率和增加了偏差,没有一个窗函数能使估计的功率谱在方差、偏差和分辨率各个方面都得到改善,因此使用窗函数只是改进估计质量的一个技巧问题,并不能从根本上解决问题。

三、 AR模型功率谱估计

1. AR模型功率谱估计简介

AR模型功率谱估计是现代谱估计中最常用的一种方法,这是因为AR模型参数的精确估计可以用解一组线性方程(Yule-Walker方程)的方法求得。其核心思想是:将信号看成是一个p阶AR过程,通过建立Yule-Walker方程求解AR模型的参数,从而得到功率谱的估计。 由于已知的仅仅是长度有限的观测数据,因此AR模型参数的求得,通常是首先通过某种算法求得自相关函数的估计值,进而求得AR模型参数的估计值。常用的几种AR模型参数提取方法有:

1) 自相关法

假定观测数据区间之外的数据为0,在均方误差意义下使得数据的前向预测误差最小。 由此估计的自相关矩阵式正定的,且具有Toeplitz性,可以用Levison-Durbin算法求解。

2) 协方差法

不作观测数据区间之外的数据为0的假设,在均方误差意义下使得数据的前向预测误差最小。由此估计的自相关矩阵式半正定的,且不具有Toeplitz性,得到的AR模型可能不稳定。

3) 修正的协方差法

不作观测数据区间之外的数据为0的假设,在均方误差意义下使得数据的前向预测误差与后向预测误差之和最小。由此估计的自相关矩阵式半正定的,且不具有Toeplitz性,得到的AR模型可能不稳定。但得到的一阶AR模型是稳定的。

4) Burg法

在约束AR模型的参数满足Levison-Durbin递归条件的前提下,在均方误差意义下使得数据的前向预测误差与后向预测误差之和最小。得到的AR模型是稳定的,但有时可能出现谱线分裂现象。

仍然用前面的仿真信号,取AR模型的阶数p=48,用上述各种AR模型参数提取方法估计的功率谱如figure2所示。

Figure2 AR模型功率谱估计的仿真结果

2. AR模型功率谱估计与经典功率谱估计性能比较

分别采用直接法和AR模型功率谱估计中的自相关法得到的上述信号的功率谱估计如figure3所示:

Figure3 AR模型功率谱估计与经典功率谱估计比较

小结:

1) 由于AR模型是一个有理分式,因而估计出的谱要比经典法的谱平滑。

2) AR谱估计的一些方法隐含着数据和自相关函数的外推,可获得高的分辨率。

3. 关于AR模型功率谱估计总结

Figure4给出了阶数对AR模型功率谱估计的影响,采用的是AR模型功率谱估计中的Burg法。

Figure4 不同阶数的AR模型功率谱估计

小结:

阶数越高,得到的谱的分辨率也越高,但方差也越大,将会产生更多的虚假谱峰;

四、 附录

本文所用来仿真的MATLAB程序如下:

% Research On Classic PSD Estimate Methods

% Author: Chen Feiqiang

% Date: 2010-12-14

%% Generate Signal

clear all; close all; clc;

randn('state',0);

Fs = 1000; % sample frequency

t = 0:1/Fs:.3;

sigma = 1; % noise variance

x = cos(2*pi*t*200) + cos(2*pi*t*210) + sigma*randn(size(t)); % generate signal:

% cosine + noise

figure;

plot(t,x), xlabel('t'), ylabel('x');

title('Signal In Time Domain');grid on;

%% Estimate PSD using periodogram method

window = rectwin(length(x)); % window function used xn = x'.*window;

index = nextpow2(length(x));

N = pow2(index);

Xw = fft(xn, N); % do N-points FFT

Pw = Xw.*conj(Xw)/N; % Calculate PSD

k = 0:floor(N/2 - 1);

figure;

plot(k*Fs/N, 10*log10(Pw(k+1)/max(Pw)));

title('Periodogram PSD Estimate');

xlabel('Frequency(Hz)');ylabel('Normalized PSD(dB)'),grid on;

hold on

%% Estimate PSD using BT method

window_a = rectwin(length(x)); % window function for data x(n)

xn = x'.*window_a;

Rx = xcorr(xn); % auto-correlation function estimate N = length(Rx);

M = floor(N/4); % the length of smooth window

% M = 100;

window_v = rectwin(M); % smooth window for BT method

RxWin = Rx(1:M).*window_v; % smooth window multiply auto-correlation function

Pw = abs(fft(RxWin, N));

k = 0:floor(N/2 - 1);

figure;

plot(k*Fs/N, 10*log10(Pw(k+1)/max(Pw)),'r');

title('BT Method PSD Estimate');

xlabel('Frequency(Hz)');ylabel('Normalized PSD(dB)'),grid on;

%% Estimate PSD using Welch method

window = 32; % the length of each segment

noverlap = 8; % overlap number for each segment nfft = pow2(nextpow2(length(x))); % nfft-points FFT for each segment

[Pxx,f] = pwelch(x,window,noverlap,nfft,Fs); % call pwelch function

figure;

plot(f,10*log10(Pxx/max(Pxx)),'g');

title('Pwelch Method PSD Estimate');

xlabel('Frequency(Hz)');ylabel('Normalized PSD(dB)'),grid on;