随机信号分析

第9章随机信号分析

随机信号和确定信号是两类性册完全不同的信号，对它们的描述、分析和处理方法也不

内向.附机信号是一种不能用确定数学关系式来描述的，无法预测未来某时刻精确值的信号，

也无法用实验的方法重复再现。

随机信号分为平稳和非平稳两类。平稳班机信号又分为各态历经和邪各态历经。本孽所

讨论的随机信号是平稳的1L是备态历经的。在研究无限长信号时,总是取某段有限长信号作

分析，这一有限长信号称为一个样本(或称子集)，而无限长信号x(D称为随机信号总体(或

林集)。各态历经的平稳随机过程中的一个样本的时间均值和集的平均值相等。因此一个样

本的统计特征代表了随机信号总体，这使得研究大大简化。］程上的随机信号一般均按各态带格品：字体颜色：红色)

历经平稳随机过程来处理。

仅在离散时间点上给出定义的随机信号称为离散时间随机信号，即随机信号序列。随机

信号序列可以是连续随机信号的采样结果，也可以是自然界里实际存在的物理现象，即它们

本身就是离散的。

平稳随机过程在时间上是无始无终的，即其能量是无限的.本身的Fourier变换也是不

存在的：但功率是有限的.j由常用功率谱密更来描述前机信号的域域特征.这是二个纭计平带格式的：字体颇色：红色—)

均的频谱特性。

平稳随机过程统计特征的计算要求信号x(n)无限长，而实际上这是不可能的，只能用

一个样本，即有限长序列来计莫.因此得到的计算也不是随机信号真正的统计值，而仅仅是

一种估计。

本章首先介绍随机信号的数字特征，旨在使大家熟悉描述甑机信号的常用特征量。然后

介绍描述信号之间关系的相关函数和协方差。这些是数字信号时间域内的描述。在频率域内，

本章介绍功率谱及其估计方法，并给出了功率谱在传递函数估计方面的应用.最后介绍描述

频率域信号之间关系的函数-一相干函数。

9.1随机信号的数字特征

9.1.1均值、均方值、方差

若连续随机信号x(t)是各态历经的.则随机信号x(t)均值可表示为：

张(3=〃*=limJC⑺力<9-i>

均值描述了随机信号的静态(直流)分量,它不随时间而变化0

随机信号X(t)的均方值表达式为：

力(9-2)

/fC/

均方位，表示了信号的强度或功率。

班机信号x(t)的均力根值衣示为：

-247―

&={l,i中(。2⑺力<9-3)

猊也是信号能量或强度的一种描述。

随机信号x（t）的方差表达式为：

E【（x-〃j］=W=啊就［小）-〃/曲（9-4）

是信号幅值相对于均值分散程度的一种表示，也是信号纯波动（交流》分属大小的

反映,

随机信号x（t）的均方差（标准差）可表示为：

a=如口北卜⑺一〃丁力（9-5）

其意义与<7；相同。

9.1.2离散随机信号

若x（n）是离散的各态历经的平稳随机信号序列。类似连续随机信号，其数字特征可由

下面式子来衣示：

均值：

年（〃）］=%=lim：£x（”）（9-6）

•V->xZa。

均方值：

£卜（”）］=《=limg次/（"）（9-7）

方差：

£■"）-凡丹=珑=］im1次卜（〃）-〃1（9&

NtbNnF

式中，内为均方根值，为均方差值。

9.1.3估计

以上计算都是对无限长信号而言，而工程实际中所取得信号是有限长的，计算中均无法

1RT―>*©或Nf+00。

对于行限长模拟随机信号，计算均值时将（9-1）式改写为：

巾"）］=凡=北必）力（9-9）

大中，均值"仅仅是对人的估计。当T足够长时，均依估计也能精确逼近真实均值出。

-w-

对于周期信号，常取T为一个周期，估计均值也就能完全代表真实均值〃L

对了有限长随机信号序列，计算其均值估计可将(9-6)式改写为：

张(")]=.="g双〃)<9-10)

当序列长度足够长时，均值估计也能精确遇近真.实均值。

类似地，可以写出均方值和方差估计表达式。

在MATLAB工具箱中,没有专门函数来计算均值、均方值和方差。但随机信号的统计数

字特征色计算都可以通过MATLAB编程实现，在数位计算中，常将连续随机信号褥散化,当

作随机序列来处理.

【例9-1】计算以长度N=100000的正态分布高斯随机信号的均值、均方值、均方根值、方

差和均方差.

%Samp9」

N=100000;%数据个数

randn('state'.O):%设置产生随机数的状态

y=randn(I,N)：%产生一个随机序列

」即('平均值:')：

yMean=mean(y)%计算防机序列的均值

dispC均平方值:‘);

y2P3*y'/N与计算其均方值,这里利用「矩阵相乘的算法

dispC均平方根:’)；

ysq=sqrt(y2p)与计综其均方根值

dispC标准差:');

ystd=std(y,+0)%相当于ystd=sqrt(sum((y-yMean).~2)/(¥T))

dispC方差:');

yd=ystd.*ystd

该程序的运行结果为：

平均值：yMean=0.0032;均平方值：y2p=1.0073

均平方根：ysq=1.0037

标准差：ystd=1.0037

方差：yd=1.0073

注意，函数s=std(x,flag)计第标准差.x为向斌或矩阵：s为标准差；flag为控制符，

用来控制标准算法。当flag=O(或缺省)时，按下式计算无偏标准差：

I

1N.3

s=-----y(x-/r)；*(9-11)

[N-yx，J

当flag=l时，按F式计算有偏标准差：

I

§=口为(演-4)叶(9-12)

9.2相关函数和协方差

9.2.1相关函数和自协方差

对于随机信号x(t),自相关函数为:

%(T)=七卜。)必+T)]=1im)，+T)力<9-13)

苴中，r为时移。

若去抑ND的均值部分，则相应的自相关函数称为自协方差，即:

Cs(?)=E{[x(r)-+T)--(9-14)

C.v(0)=a；(9-15)

对于离散随机信号序列，x(n)的自相关函数和自协方差为：

IA'-l

(("[)=E[.v(n).r(ZJ-w)]=limTvX+〃])(9-16)

Nf*Nn.Q

.(〃?)=七心(〃)一从卜("+〃】)—〃』｝=/?«〃)—〃：<9-17>

式中，m为延迟。

9.2.2互相关函数和互协方差

对于两个不同随机信号x(t),y(t),互相关函数为:

&,(，)=4v")N，+，；]=lim工)力<9-18)

式中，r为时移。

互协方差为：

对于互协方差有：

&,(r)=&O(9-20)

对于离散班机序列x(n)和y(n),互相关函数和互协方差为：

&、.(〃?)=E[x(〃).y("+⑼]=1im：£.«")》(〃+〃i)(9-21)

Nf®Nn=0

&(，〃)=E&〃)-〃IN”+切)-〃、])=勺(2(9-22)

［带格式的：字体颜色：红色

注我到上面的公式大要求.rfoo/XN->8j而在匚程中信号是有限长的，因此只能［带格式的：字体颜色：红色1

得到相关函数和协方至的估计值，当t或N足够长时，估计值能精确地遍近其实值。「带格式的:『体就色：红色:

〔用格式的：字体颇色：红也】

带格式的：字体颜色：红色:

-247-

9.2.3计算相关函数和协方差的MATLAB函数

VATLAB信号处理.1：具箱提供「计算班机信号相关函数xcorr。

函数xcorr用于计算随机序列自相关相互相关函数。调用格式为：

[c,lags]=xcorr(x.y[,maxlags,,option'])

式中，x.y为两个独立的随机信号序列，长度均为N：c为x,y的互相关估计：lags为相关估

计c的序号向量，其范围为[-maxlagscmaxlags]。

oplion缺省或‘none’时，函数xcorr按下式执行非归一化的相关：

N-H-1

"此°

n=o(9-23)

c*(-WI)rn<0

上角加星号的变最衣示其共施。

option为‘biased’,计算有偏互相关函数估计

=="，〃)(9-24)

option为'unbiased',计算无偏互相关函数估计i至格式的：字体献色：红色一

C^(m)=-^CKV(m)(9-25)

N-网

option为'coeff',序列归一化,使零延迟的自相关函数为1。

maxlags为x和y的最大延迟，该项缺省时，函数返回(ftc的长度为2NT：该项不㈱省时，

函数返【可值c的长度为2*maxlngs+1.

该函数也可用于求一个随机信号序列x(n)的自相关函数，调用格式为：

c=xcorr(x,maxlags)

MATLAB信号处理工理箱还提供了计免协方差方数XCOY.若x为■个向艮，调用LXMY(X)

函数，则输出c为x序列的方差：若x为一个矩阵，则c=xcov("输出更阵x的协力差起阵,

其对角线元,素为每一列数据的方差,芈对角线元素.表东对应行和对应列的协方弟。调用

c=xcov(x,y),其中凡丫具仃相同的长度,则c为式(9-22)所对应的互协方差。

【例9-2]求带有白噪声干扰的频率为101k的正弦信号和白噪声信号的自相关函数为进行

匕较。

%Samp9_2

elf;N=1000：Fs=500：家数据长度和采样频率

n=0:N-l：t=n/Fs;%时间序列

Lag=100;舟延迟样点数

randn('state',0);%设置产生随机数的初始状态

x=sin(2*pi»10*t)+0.6*randn(l,length(t))：%原始信号

[c,lags]=xcorr(x,kig,,unbiased')：'对原始信号进行无偏自相关估计

subplot⑵2,1),plot(t,x);%绘原始信号x

xlabelC时间/s');ylabel('x(t)");ti:le(,带噪声周期信号');gridon：

subplot(2.2,2)：plot(lags/Fs,c)：'绘x信号自相关Jags/Fs为时间序列

xlabel('时间/s');ylabelCRx(t),);

title('带噪声周期信号的自相关')；gridon;

、信号xl

xl=randn(Llength(x));为产生一与x长度一致的随机信号

[c,lags]=xcorr(xl.Lag,*unbiased');%求随机信号xl的无偏自相关

subplot(2,2,3),plot(t.xl)：%绘制随机信号xl

xlabelC时间/s');ylabel('xl(t)');title('噪声信号');

gridon;

subplot(2.2,4)：plot(lags/Fs,c)：为绘制随机信号xl的无偏自相关

xlabelC时间/s')：ylabel('Kxl(t),):

title('噪声信号的自相关');

gridon

带噪声周期信号的自相关

时间/S

图％I例9-2求得带有白曝尹的周期信号和噪声信号的自相美

程序运行结果见图9-1。可以看出，含有周期成分和噪声成分的自相关函数在7=0时

具有最大值。且在r较大时仍具有明显的周期性，其频率和周期信号的频率相同；而不含周

期成分的纯噪声信号仕：■=()时具仃最大值.且在r稍大时很快衰减至冬。自相关的这一性

J4T

质被用来识别随机信号中是否含有周期成分并用于确定所含周期成分的频率。

【例9-3】已知两个周期信号：x(r)=sin(2磔),y(r)=O.5sin(2磔+90"),其中f=：0Hz,

求互相关函数兄,卜)。

%Samp93

Cir：x=iooo；rs=500；、数据长度和采样频率

n=0:N-l;t=n/Fs：%时间序列

Lag=200：%最大延迟单位

x=sin(2*pi*10*t)：气周期信号x

y=0.5*sin(2*pi*10*t-90*pi/180)：%与x有90'相移的信号y

[c,lags]=xcorr(x,y.Lag,'unbiased');%求无偏互相关

subplot.⑵1,1):plot(t,x,'r)：―绘制x信号

holdon：plot(t,y,':')：$在同一媒图中绘y信号

legendCx信号'，'y信号')

xlabelC时间/s')jylabel('x(t)y(l)');

titleC原始信号’)；gridon：

holdoff

cnhplot(2,I,2),pInt(Iagc/l'c,<•/r')：(绘制x,y的互相关

xlabelC时间/$');ylabel('Rxy(t),);

title('信号x和y的相关'):gridon

^47-

信号X和y的相关

图9-2例％3所给的炭华信号及其互相关

程序运行结果为图9-2。由结果可见，&,(「)也是周期为10Hz的余弦信号(注意两个

图形其中的时间尺度有改变)，其幅位为0.25,初始相位为90'。

此可得到“相关函数的一个虫要特性：两个均值为零IU!行相同频率的周期信：，儿1带格式的：字体颜色：红色

互相关函数R“(r),保留原信号频率、相位差和幅内信息.:嚣鬻蓼器嚣

下面例子给出互相关函数的另一个工程应用实例“若信号x(t)和信号y(i)是同类信号

且有时移，用互相关函数可以准确地计算出两个信号时移大小。这种特性使得互相关函数广

泛地应用于工程测量技术中。

【例9-4】两个sine信号有0.2s的时移,用互相关函数计算时移的大小。

%Siunp9_4

elf

N=1000:n=0:N-l:Fs=500:t=n/Fs;%数据个数采样频率和时间序列

Lag=200;%最大延迟单位数

xl=90*sinc(pi*(n-0.l*Fs));%第一个原始信号，延迟0.1s

y1=50*sinc(pi♦(n-0.3*Fs));%第二个原始信号，延迟0.3s

[c,lags]=xcorr(xl,yl,Lag,'unbiased'):'计算两个函数的互相关

subplot(2.1,l),plot(t,xl,'r);%绘第一个信号

holdon;plot(t,yl,'b/)：%在同一幅图中绘第二个信号

logendC信号x','信号y');%绘制图例

xlabelfBj(BJ/s*)：ylabel('x(t):

title('信号x和y");holdoff

subplot(2,1.2),plot(lags/Fs,c,'r'):%绘制互相关信号

xlabelC时间/s');ylabcl('Rxy(t),);

titleC信号x和y的相关')；

程序运行结果见图93。可以清楚地看到第二个信号相对于第一个信号延迟了0.2s,即

石-0.2s处出现相关极大值。因此可以采用该项技术检测延迟信号。

信号x和y

100

信号x

80信号y・

60

!，。

20

0

0.81

时间/s

信号，和y的相关

时间“

图9-3例9-4两个sine信号的互相关

J4T

9.3功率谱估计

M率谱估计的目的是根据仃限数据给？储点的机过程的顺率成分分布的描述。上节所(带格式的：字体颜色：红色

止的相关分忻是时域内在噪声背景下提取有用俗息的途径，而功率谱是姚域内提取淹没在噪

声中有用信息的分析方法。

9.3.1功率谱,密度

假如随机信号X⑴的自相关函数为RJr),R,(T)的Fourier变换为：

S,(/)=「R'(W"r<9-26)

则定义S,⑴为x⑴的自功率谱密度或称为自功率谱.尸1为S,⑹可解杼为x⑴俏平均1带格式的：字体颜色：红色)

功率相对卜频率的分布函数。自功率谱S，C)包含R,(r)的全部信息。如果噪声信号中含有

某种频率成分，可以从自功率谱中看出。

若随机信号x(i)的自功率谱为£(「).则根据Fcuriar逆变换可得•

R、(r)=[S*(/)/■“<9-27>

和自功率谱类似，两个陵机信号x(t)木1y(t)互相关的频率特性可用互功率谱密度来描

达。互功率谱密度和互相关函数也是一Fourier变换对。

S”(/)=<9-28)

R”(7)=[S.v(/)e'W(9-29)

对于离散随机序列x(n),自功率谱密度S,(f)和自相关函数R,Gn)的关系为：

X-[带格式的：字体颓色：红色____I

5;(/)=(9-30)〔带格式的：字体颜色：红色

AJ________212^_____________a___________________________1带格式的:字体颜色：红色J

这里，Ts为数据采样间隔。

对于离散随机序列x(n)和y(n),互功率谱密度S..(f)和互相关函数Rl((m)关系为：

5")=£勺(Me"—(9-31)

ms-®

旦有

2

S-D=SQ,Sx(f)Sy(f)>|Sty(/)|(9-32)

浜际工程中随机序列长度均为有限长，因此利用有限长瓶机序列计算的自功率诺密度和〔带格式的：字体颜色：红色1

互功率谱密度只是真实值的一种估计。

功率谱估计方法一般可分为伊数估计和非参数估计两类。MATLAB信号处理工具箱介绍[带格式的：字体颜色：红色)

的功率谱态度非参数估计方法有WELCH法、MTM法和MUSIC法；参数估计方法有MEM法，卜

面先弁绍求取功率谱估计的基小力法，然后逐•介绍信号处理工具箱中绐出的上述力去。

-249―

9.3.2周期图法

1.基本方法

周期图法是直接将信号的采样数据x(n：，进行/ourier变换求取功率谱密度估计的方法。

饭定有限长的机信号序列为xGD.它的Fourier变换和功率谱密度估计4(/)存在下面的

关系：

Sl(f)=jj\x(f)\'(9-33)

式中，N为随机信号序列x(n)的长度.在段散的频率点f=W，有：

1(A)=-l|X(初,=—|FF7(x(n)］\k=OA,--,N-\(9-34)

其中，l：FT［x(n)］为对序列x(n)的Fourier变换，由于FFT〔x®］的周期为N,求得的

功率谱估计是以N为周期，因此这种方法称为周期图法。下面用例子说明如何采用这种方法

进行功率谱估计。

【例9-5］用Fourier变换算法求信号x")=sin(2城/)+2sin(2M/)+的功率谱。

其中fi=50Hz,f2=120Hz,w(t)为白噪声,采样频率Fs=1000Hz0(1)信号长度N=256：(2)信

号长度N=1024。

%Samp95

clf;I's=1000;

%第一种情况：N=256

N=256:Nfft=256:%数据长度和FFT所用的数据长度

n=0:N-l：t=n/Fs：为采用的时间序列

xn=sin(2*pi*50*t)+2*sin(2*pi*120*l：-+randn(l.N);气带有噪声的信号

Pxx=10*logl0(abs(fft(xn,Nfft).2)/K):%Fourier振幅谱平方的平均值,并转换为dB

f=(0:1ength(Pxx)-1)*Fs/length(Pxx);%给出频率序列

subplot(2,1,1),plot(f,Pxx);%绘制功率谱曲线

xlabelC频率/Hz');ylabel(,功率谱/dB')：

title('周期图N=256');gridon;

为第二种情况：N=1024

Fs=100();N=1024：Nfft=1024;%采样频率、数据长度和IFF所用数据长度

n=0:N-l;t=n/Fs;*时间序列

xn=sin(2*pi*50*t)+2*sin(2*pi*120*t?+randn(1,N);舟带有噪声原始信号

Pxx=10*logl0(abs(fft(xn,Nfft).*2)/N):%Fourier振幅谱平方的平均值,并转换为dB

f=(0:1ength(Pxx)-1)*Fs/1ength(Pxx)：、频率序列

subplot(2,1,2),plot(f,Pxx);%绘制功率谱曲线

xlabelC频率/Hz');ylabel('功率谱/dB');

title、周期图N=1024')：gridon

程序的运行结果为图9-4,可以看出，右频率50Hz和120Hz处，功率谱有两个峰值，说

中信号中含有50Hz和120Hz两种频率成分,但功率谱密度在很大范围内波动，而且广没有

国信号取样点数N的增加加有明显改进。注意本例中500Hz为Nyquist频率，由第3章介绍

^4^

ftjFourier变换的知识可知，只观看500Hz之前的功率谱可以推知以后的功率谱。在本隼后

面的例题中我们将只研究Nyquist频率之前的功率谱。

8

5

生

阴

日

图94例9-5中含有噪声信号的功率诺

用有限长样本序列的Fourier变换来表示随机序列的功率谱，只是一种估计或近似，不

可避免存在误差。力/她少法入使功的含侑计史如平滑，”J采用分段平均周期图法(带格式的：字体颜色：红色

(Bartlett法人加窗平均面期图法(Welch法)等方法加以改进。

2.分段平均周期图法(Bartlett法)

将信号序列x5),n=0,分成互不重序的P个小段,每小段由m个采样值，

则1加=%对每个小段信号序列进行功率谱估计,然后再取平均作为整个序列x(n)的功率谱

估计.

平均周期图法还可以对信号x(n)进行亘段分段，如按2:1重总分段，即前一段信号和

后•段信号有•半是重登的。对母•小段信号序列进行功率诺估计，然后再取平均值作为整

-247-

个序列x(n)的功率谱估计。

这两种方法都称为平均周期图法，•般后者比前看好°

【例9-6】对例9-5中的信号序列，用两种平均周期图法求信号的功率谱密度估计。

%Samp96

clf;Fs=1000;%数据采样点数

%运用信号不重登分段估计功率谱

N=1024：Nsec=256;n=0:N-1：i=n/Fs：%数据点数.分段间RS，时间序列

randnCstale',0):先设更产生随机数的初始状态

xn=sin(2*pi*50*t)+2*sin(2*pi*120*t?+randn(LN);%带噪声的原始信号

pxxl=abs(fft(xn(l:256),Nsec).'2)/Nsec：%第一段功率诺

pxx2=abs(fft(xn(257:512),Nsec).*2)^sec;%第二段功率谱

pxx3=abs(fft(xn(515:768),Nsec).*2),/Nsec;$第三段功率谱

pxx4=abs(fft(xn(769:1024),Nsec).*2j/Nsec;%第四段功率谱

Pxx=10*1og10((pxx1+pxx2+pxx3+pxx4)/4)：气平均得到整个序列功率谱

t=(0:1ength(l>xx)-1)»I's/length(Pxx)：益给出功率谱对应的频率

subplot(2,1,l)tplot(f(l:Nsec/2),Pxx(1:Nsec/2))：舟绘制功率谱曲线

xlabelC频率/Hz')：ylabel('功率谱/iff);

titleC平均周期图(无重叠)N=4*256，)：

gridon

%运用信号重登分段估计功率谱

pxxl=abs(fft(xn(1:256),Nsec).-2)/Nsec：机第•段功率诺

pxx2=abs(fft(xn(129:384),Nsec)."2)fUsec;先第二段功率谱

pxx3=abs(ffI(xn(257:512),Nsec).*2)/Nsec;%第三段功率谱

pxx4=abs(fft(xn(385:640),Nsec).2)/Nsec;$第四段功率谱

pxx5=abs(fft(xn(513:768),Nsec).-2)/Nsec：%第五段功率谱

pxx6=abs(fft(xn(641:896),Nsec).-2)/Nsec：先第六段功率谱

pxx7=abs(ffI(xn(769:1024),Nsec).*2)/Nsec;%第七段功率谱

Pxx=10*logl0((pxxl+pxx2+pxx3+pxx4*pxx5+pxx6+pxx7)/7)：%功率谱平均并转化为dB

f=(0:length(Pxx)-l)*Fs/length(Pxx):%频率序列

subplot(2,1,2),p1ot(f(1:Nsec/2),Pxx(1:Nsec/2))：%绘制功率谱曲线

xlabelC频率/Hz');ylabelC功率谱/dB');

title{'平均周期图(重看一半)N=1024')：

gridon

平均周期图（无亚0）%4256

平均周期图（重®一半）N=1024

电

、

躯

脩

屋

图9-5例9-6中重叠和木页无分段计匏的功率谐比较

程序运行结果为图95,上图采用不重注分段法的功率谱估计，下图为2:1重登分段的

功率谱估计，可见后者估计曲线较为平滑。与上例比较，平均周期图法功率谱估计具有明显

效果（涨落曲线靠近OdB）。

3.加窗平均周期图法

加窗平均周期图法是对分段平均周期母法的改进。在信号序列x（n）分段后，用非矩形

窗口对每一小段信号序列进行预处理，再采用前述分段平均周期图法进行整个信号序列x（n）

的功率谱估计。由窗函数的基本知识（第7率）可知，采用合适的非矩形窗口对信号进行处

理可减小“频漕泄露”，同时可增加频峰的宽度，从而提高频谱分辨率。

【例9-7】对例9-5中的信号序列.用加腐平均周期图法进行功率谱密度估计。

%Samp9_7

clf;Fs=1000;%采样频率

N=1024;Nsec=256：n=0:N-l;t=n/Fs;为数据长度,分段数据长度、时间序列

w=hanning（256）';%采用的窗U数据

randnCstate,,0);%设置产生随机数的初始状态

xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N)：%带噪声的信号

%采用不重叠加窗方法的功率谱估计

pxxl=abs(fft(w.*xn(l:256),Nsec).*2)/norm(w)2：%第一段加窗振幅谱平方

pxx2=abs(fft(w.*xn(257:512),Nsec).-2)/norm(w)'2;%第二段加窗振幅谱平方

pxx3=abs(fft(w.*xn(513:768),Nsec).2)/norm(w)'2：%第三段加窗振幅谱平方

pxx4=abs(rn(w.*xn(769:1024)，Nsec)/2)/norm(w>2;%第四段加窗振幅谱平方

Pxx=10*logl0((pxxlfpxx2+pxx3*pxx4)/4)：先求得平均功率谱，转换为dB

f=(0:1ength(Pxx)-1)*Fs/length(Pxx):%求得频率序列

subplot(2,1,1),plot(f(1:Nsec/2),Pxx(l:Nsec/2))：%绘制功率谱曲线

xlabelC频率/Hz');ylabel('功率谱/dB');

ti11eC加窗平均周期图(无电桂)N=4*256J);

gridon

外采用重登加窗方法的功率谐估计

pxxl=abs(fft(w.*xn(1:256),Nsec).*2)/norm(w)'2：%第一段加窗振幅谱平方

pxx2=abs*xn(129:384),Nsec).*2)/nor»(w)*2;%第二段加窗振幅谱平,方

pxx3=abs(fft(w.*xn(257:512),Nsec).2)/norm(w)'2;*第三段加窗振幅谱平方

pxx4=abs(fft(w.*xn(385:640),Nsec).2)/norm(w)-2;为第四段加畲振幅谱平方

nxx5=«bs(fft(»-*xn(513t768),Nsec)2)/nora(w)*2：%第五段加窗振幅谱平方

pxx6=abs(fft(w.*xn(6-11:896),Nsec).2)/norm(w)*2;舟第六段加窗振幅谱平,方

pxx7=abs(fft(w.*xn(769:1024),Nsec).'2)/nonn(w)'21%第七段加窗振幅谱平方

Pxx=10*1og10<(pxx1+pxx2+pxx3+pxx4+pxx5+pxx6+pxx7)/7):%平均功率谱转换为dB

f=(0:length(Pxx)-1)*Fs/length(Pxx);%频率序列

subplot(2,1,2),plot(f(l:Nsec/2),Pxx(1:Nsec/2));、绘制功率谱曲线

xlabelC频率/Hz');ylabcl('功率谱/dB')：

title/加窗平均周期图(重叠一半)N=1024');

gridon

程序运行结果为图9-6,其中上图采用无重登数据分段的加窗平均周期图法进行切率谱

估计，而下图采用重叠数据分段的加窗平均周期图法进行功率谱估计，显然后者是更佳的.

信号漕峰加宽,而噪声谱均在OdB附近，更为平坦（注意采用无重叠数据分段噪声侑最大

的下降分贝数大于5dB,而亚比数据分段周期图法噪声的最大下降分贝数小于5dB）,

^47-

加密平均周期图（无*0）N=4,256

加窗平均周期图便箴一半）N=1024

物

协

£

图9-6例9-7给出的不正登加密和重叠加两周期图法得到的功率谱曲线

「格式的：④色：红色

4.Welch法估计及其MATLAB函数

Welch功率谱密度就是用改进的平均周期图法来求取随机信号的功率谱密度估计的.

Welch法采用信号1R登分段、加窗函数和l'1'r算法等”•第一个信号序列的自功率谱估i（PSD

如上例中的下半部分的求法）和两个信号序列的互功率谱估计（CSD）,

MATLAB信号处理工具箱函数提供了专门的函数PSD和CSD自动实现肥kh法估十，而

不需要像以上例子一样自己编程。

（1）函数psd利用Welch法估计一个信号自功率谱密度,函数调用格式为:

[Pxx[,f]]=psd（x[,Nfft,Fs,window,Noverlap,'dflag*]）

式中,x为信号序列：Nfft为采用的FFT长度。这一值决定了功率谱估计速度,当Nfft采

用2的耗时,程序采用快速免法：心为采样领率；Wind。”定义窗函数和x分段序列的长度。

窗函数长度必须小于或等于Nfft,否则会给出错误信息："verlap为分段序列于叠侑采样

点数（长度），它应小于Nfft;dfla*为去除信号趋势分量的选择项：'linear',去除线性

电势分量mean'去除均值分量，‘none'不做去除趋势处理。Pxx为信号x的自勺率谱

-247-

密度估计。f为返回的频率向量，它和Pxx对应，并且有相同长度。

在psd函数调用格式中，缺省值为：Nfft=min(256,length(x)),Fs=2Hz,

window=hanning(Nfft),noverlap=O.若x是实序列，函数psd仅计算频率为正的功率谱估

计。

【例9-8】对例9-5中的信号序列，用函数psd绘制自功率谱密度估计曲线，

%Samp9_8

cir：Fs=1000：%采样领率

N=1024;Nfft=256;n=0:N-l;t=n/Fs;MS据长度、时间序列

window=hanning(256);，选用的窗口

noverlap=128;%分段序列重登的采样点数(长度)

dflag=，none';%不做趋势处理

randn(，state',0);%设置产生随机教的初始状态

xn=sin(2*pi*50*t)+2*sin(2*pi*120*t?+randn(1,N)：与带噪声原始信号

Pxx=psd(xn,Nfft,Fs,window,novorlap,dflag)：%功率谱估计

f=(0:Nfft/2)»l-s/Nfft：%求得对应的频率向量

plot(f,10*logl0(Pxx))：%绘制功率谙

xlabelC频率/Hz');ylabelC功率谱/dB');

titleCPSD—Welch方法')：gridon

程序运行结里见图9-7.可以看到.采用Welch法求得的功率谱与前面所述方法相比效

果最好。使信号得以突出，而抑制了噪声。

J4T

PSD-Welch方法

25

_10"，________：_______：______：______：______!______：_______：______•

050100150200250300350400450500

频率，'Hz

图9-7Welch法得到的功率诺估计

注意程疗前半部分中频率向量f的创建方法。它与函数psd的输出Pxx长度的关系如下:

若x为实序列,当Nfft为奇数时,f=(0:(Nfft+l)/2-l)/Nfft：SNfft为偶数时，

f=(0:Nfft/2)/Nfft.

函数还有一种缺省返回值的调用格式，用于直接绘制信号序列x的功率谱估”•曲浅。

函数还可以计算带有置信区间的功率谙估计，调用格式为：

[Pxx,Pxxc,f]=psd(x,Nfft,Fs,window,Noverlap,p)

式中，p为置信区间，

【例9-9】设计一个归一化频率为0.2的FIR滤波洛，对一个白噪声信号序列进行滤波，对

滤波后的信号绘制置信区间为0.95的功率谙估计曲线。

%Samp9_9

clf;h=firl(30,0.2,boxcar(31)):%设计30阶截止频率为0.2的FIR1逑波器

r=randn(16384,1);“产生白噪声随机向量

x=filtor(h.1,r)：%对随机向发泄波

^47-

psd(x,1024,10000,kaiser(512,5),0,0.95);%采用Kaiser窗进行功率诺估计

xlabeK'频率/Hz)：ylabel(：功率谱SB')

捏序运行结果见图9-8.程序中首先产生一个FIR谑波器，其梭止频率为0.2,和对于

本例中10000Hz的采样频率，5000Hz对应的妇一化频率为1.0.2的归一化频率对应于

1300Hz。噪声为白噪声信号，含有多种频率成分。因此采用上述滤波涔对噪声进行滤波后,

截止频率1000Hz以后的信号大大衰减，具“丰富频率成分的白噪声经过滤波器泄波后输出

信号领增相当于渔波器的幅坝特性曲线。其中上面曲战和下面的纹分别对成丁一95%的置信区

间的上下限。

B

P

、

物O

-3

SHS-

-70

0500100015002000250030003500400045005000

频率/Hz

图9-8例9-9设计注波器对口噪，”泄波后的功率谐

此可知.•沙波器输入臼噪声序列的输出伉漕或门相关可以确定滤波器俏拨率[带格式的：字体颜色：红色

特性.

(2)函数csd利用welch法估计两个信号的互功率谱密度，函数调用格式为:

[Pxy[,f]]=csd(x,y,Nfft,Fs,window,Noverlap,fdflag')

[Pxy,Pxyc[,f]]=csd(x,y,Nfft,Fs,window,Noverlap,p)

-247-

这里，x,y为两个信号序列：Pxy为x,y的W功率谱估计；其他参数的意义同自功率谙函数

psdo

【例9-10】产生两个长度为16384的白噪声信号和一个带有白噪声的1000Hz的周期信号,

求两个白噪声信号及白噪声与带有噪声的周期信号的互相关谱，并绘制互相关谱曲线，FFT

所采用的长度为1024,采用500个点的三角形窗，并且没有重登，采样频率Fs=10000Hz。

%Samp910

Fs=10000：%信号采样领率

randnCslate',0):*设置随机状态初始值

x=randn(16384,1);%产生第一个白噪声信号

y=randn(16384,1)：舟产生第二个白噪声信号

z-2*sin(2*pi*1000*[1:16：i84]*/Fs)+ran<ln(l6384,1):%产生含周期成分的噪声信号

[Pxy,f]=csd(x,y,1024,Fs,triang(50Q),0);%白噪声信号互相关谱，无束叠三角窗

[Pxz,f]=csd(x,z,1024,Fs,triang(50C),0)个噪声与周期信号互相关谱无重叠三角窗

subplot(2,1,l),plot(f,10*logKXPxy)),gridon;先绘制x,y信号互相关谱

ylabel('互相关谱/dB'),title('噪声信号的互相关谱')

subplot(2,1,2),plot(f,10*logl0(Pxz)),gridon;%绘制x,z信号的互相关谱

ylabelC互相关谱/dB'),xlabel('频率/Hz')

titleC带噪声周期信号的互相关谱')

程序运行结果为图9-9.可以看到.两个白唾声信号的互功率谱(上图)杂乱无序.若

不出周期成分，大部分功率谱在-5dB以卜然而白噪声与带有噪声的周期信号的功率谱在

其周期(频率为1000Hz)处有峰值，清焚地表明了周期信号的周期或频率。因此，利用

未知信号与白噪声信号的互功率谱也可以检测未知信号中所含有的频率成分.

J4T

0具声信号的互相关谱

CPD

、

期

*

2

101

CPD

、

*轴

耍

Ini

图99例9-10的两个噪声信