功率谱分析

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

姓名:***

学号：************

专业：电子信息工程

一．原理

古典谱估计之相关函数法

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

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

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

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

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

一次是估计时仅利用了从-(M-1)到(M-1)的,这相当于对加矩形窗截断,将截成(2M-1)长,这称为加延时窗.式中和分别表示对和的估值.由于M ?N,使得上式的运算量不是很大,因此在FFT

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

古典谱估计之周期图法

首先，数据加窗得到有限长序列x N (n)直接求傅里叶变换，得频谱X N (e i ω),即

X N e iω = x N (n )e

?iωn

N?1n =0 然后，取频谱幅度的平方，并除以N ，以此作为对x(n)真实功率谱Sx(e i ω)的估计，即

Sx e i ω =1

|X N e iω |2

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

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

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

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

N 1

x N N n 0

1R (m)x (n)x (n m)N -==+∑1(1)

??S ()()M j j m x

x

m M e R m e ωω--=--=∑

x R (m)()N x n ?S ()j x e ωx R (m)x R ()m x R ()m x R (m)?S ()j x

e ωx R ()m S ()j x e ω

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

可以看出，Ｐ阶ＡＲ模型有Ｐ＋１个待定系数，由上式，可得白噪声激励得到的系统输出

可以理解为，用ｎ时刻之前的ｐ个值的线性组合

来预测ｎ时刻

的值预测误差为．在均方误差最小准则下，组合系数

的选

择应使预测误差的均方值最小．经过一系列的运算，最终可以得到ＡＲ模型的正则方程

也就是尤勒－沃克（Yule-Walker ）方程．

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

再由公式Sy e iω =Sx e iω ?|H (e iω)|2可以得出功率谱。

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

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

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

可以看出，Ｐ阶ＡＲ模型有Ｐ＋１个待定系数，由上式，可得白噪声激励得到的系统输出

1G ()1m

i

i i H z a z -==

+∑1

()()()

p

i i x n a x n i G n ω==--+∑1

()

p

i i a x n i =--∑()x n ()G n ω12,,,p

a a a ()G n ω1

21(),1,2,()(),0p

i x i x p

i x i a R m i m p R m a R i G m ==?--=??=??-+=??∑∑1G ()1m

i

i i H z a z -==

+∑1

()()()

p

i i x n a x n i G n ω==--+∑

可以理解为，用ｎ时刻之前的ｐ个值的线性组合

来预测ｎ时刻

的值预测误差为．在均方误差最小准则下，组合系数

的选

择应使预测误差的均方值最小．经过一系列的运算，最终可以得到ＡＲ模型的正则方程

也就是尤勒－沃克（ｙｕｌｅ－ｗａｌｋｅｒ）方程．

在上面方程中令ｐ＝１，得一阶ＡＲ模型的尤勒-沃克方程为

可以解得

进而求得

在一阶的基础上进行递推，将阶次为ｍ时的第ｍ个系数定义为反射系数，用表示，于是可以将计算ｍ阶ＡＲ模型参数的Levinson-durbin 递推算法表示如下：

式中，．

1

()

p

i i a x n i =--∑()x n ()G n ω12,,,p

a a a ()G n ω1

21(),1,2,()(),0p

i x i x p

i x i a R m i m p R m a R i G m ==?--=??=??-+=??∑∑111(0)(1)(1)(1)(1)(0)0x x x x R a R R a R ρ+=??

+=?1(1)(1)(0)

x x R a R =-

2211(1)

(0)(0)[1(1)]

(0)x x x x R R R a R ρ=-=-()m a m m k 1

11

1112

1[()()()]()()(),1,2,3,,1(1)m x m x i m m m m m m m m m R m a i R m i k a i a i k a m i i m k ρρρ--=----?

-+-?

=????

=+-=-???

?=-???∑1,2,

m p =2

0[()](0)x E x n R ρ==

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

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

ε=12 εf +εb =12

[e m?1f n ]2N?1

n =m + [e m?1b n ?1 ]2N?1

n =m

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

e m f

n =e m?1f

n +k m e m?1b

n ?1 e m b n =e m?1

b n ?1 +k m e m?1f n e 0f

n =e 0b n =x (n ) 得到反射系数的估计公式为：

k m =2 [e m?1f

n e m?1fb

n ?1 ]

N?1

n =m [e m?1

f n ]2+[e m?1

b n ?1 ]2 N?1n =m

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

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

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

Levinson 关系式如下：

1

-m 1,2,i i),-(m (i)(i)1-m 1-m m

，?=+=a k a a

m

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

ρ

m

2

G =

其中

ρ

m

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

)

-(12

m 1

-m m

k ρ

ρ

= 求得。

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

ρ

ρ0

=[

]

(0)

(n)E x 2

R x =

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

)(e j x ω

S =2

j 2)H(e ωωσ

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

二．程序，图像及分析

时域波形图

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(e iω)的一致估计。于是我们利用平滑和平均的方法改进周期图法。

图三

图四

图五

图六

平滑是用一个适当的窗函数W(e iω)与算出的功率谱进行卷积，使谱线平滑。这种方法得出的谱估计是无偏的，方差也小但分辨率下降。平均是将截取的数据段x N(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+p

N?p ρpρp随阶次增加而减小，而N+p

N?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）