9随机信号分析.docx
《9随机信号分析.docx》由会员分享，可在线阅读，更多相关《9随机信号分析.docx（45页珍藏版）》请在冰豆网上搜索。
9随机信号分析
第9章随机信号分析
随机信号和确定信号是两类性质完全不同的信号，对它们的描述、分析和处理方法也不相同。
随机信号是一种不能用确定数学关系式来描述的，无法预测未来某时刻精确值的信号，也无法用实验的方法重复再现。
随机信号分为平稳和非平稳两类。
平稳随机信号又分为各态历经和非各态历经。
本章所讨论的随机信号是平稳的且是各态历经的。
在研究无限长信号时，总是取某段有限长信号作分析。
这一有限长信号称为一个样本（或称子集），而无限长信号x（t）称为随机信号总体（或称集）。
各态历经的平稳随机过程中的一个样本的时间均值和集的平均值相等。
因此一个样本的统计特征代表了随机信号总体，这使得研究大大简化。
工程上的随机信号一般均按各态历经平稳随机过程来处理。
仅在离散时间点上给出定义的随机信号称为离散时间随机信号，即随机信号序列。
随机信号序列可以是连续随机信号的采样结果，也可以是自然界里实际存在的物理现象，即它们本身就是离散的。
平稳随机过程在时间上是无始无终的，即其能量是无限的，本身的Fourier变换也是不存在的；但功率是有限的。
通常用功率谱密度来描述随机信号的频域特征，这是一个统计平均的频谱特性。
平稳随机过程统计特征的计算要求信号x（n）无限长，而实际上这是不可能的，只能用一个样本，即有限长序列来计算。
因此得到的计算值不是随机信号真正的统计值，而仅仅是一种估计。
本章首先介绍随机信号的数字特征，旨在使大家熟悉描述随机信号的常用特征量。
然后介绍描述信号之间关系的相关函数和协方差。
这些是数字信号时间域内的描述。
在频率域内，本章介绍功率谱及其估计方法，并给出了功率谱在传递函数估计方面的应用。
最后介绍描述频率域信号之间关系的函数---相干函数。
9.1随机信号的数字特征
9.1.1均值、均方值、方差
若连续随机信号x（t）是各态历经的，则随机信号x（t）均值可表示为：
（9-1）
均值描述了随机信号的静态（直流）分量，它不随时间而变化。
随机信号x（t）的均方值表达式为：
（9-2）
均方值
表示了信号的强度或功率。
随机信号x（t）的均方根值表示为：
（9-3）
也是信号能量或强度的一种描述。
随机信号x（t）的方差表达式为：
（9-4）
是信号幅值相对于均值分散程度的一种表示，也是信号纯波动（交流）分量大小的反映。
随机信号x（t）的均方差（标准差）可表示为：
（9-5）
其意义与
相同。
9.1.2离散随机信号
若x（n）是离散的各态历经的平稳随机信号序列。
类似连续随机信号，其数字特征可由下面式子来表示：
均值：
（9-6）
均方值：
（9-7）
方差：
（9-8）
式中，
为均方根值，
为均方差值。
9.1.3估计
以上计算都是对无限长信号而言，而工程实际中所取得信号是有限长的，计算中均无法取
或
。
对于有限长模拟随机信号，计算均值时将（9-1）式改写为：
（9-9）
式中，均值
仅仅是对
的估计。
当T足够长时，均值估计
能精确逼近真实均值
。
对于周期信号，常取T为一个周期，估计均值
就能完全代表真实均值
。
对于有限长随机信号序列，计算其均值估计可将（9-6）式改写为：
（9-10）
当序列长度足够长时，均值估计
能精确逼近真实均值
。
类似地，可以写出均方值和方差估计表达式。
在MATLAB工具箱中，没有专门函数来计算均值、均方值和方差。
但随机信号的统计数字特征值计算都可以通过MATLAB编程实现。
在数值计算中，常将连续随机信号离散化，当作随机序列来处理。
【例9-1】计算以长度N=100000的正态分布高斯随机信号的均值、均方值、均方根值、方差和均方差。
%Samp9_1
N=100000;%数据个数
randn（'state',0）;%设置产生随机数的状态
y=randn（1,N）;%产生一个随机序列
disp（'平均值:
'）;
yMean=mean（y）%计算随机序列的均值
disp（'均方值:
'）;
y2p=y*y'/N%计算其均方值,这里利用了矩阵相乘的算法
disp（'均方根:
'）;
ysq=sqrt（y2p）%计算其均方根值
disp（'标准差:
'）;
ystd=std（y,1）%相当于ystd=sqrt（sum（（y-yMean）.^2）/（N-1））
disp（'方差:
'）;
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=0（或缺省）时，按下式计算无偏标准差：
（9-11）
当flag=1时，按下式计算有偏标准差：
（9-12）
9.2相关函数和协方差
9.2.1相关函数和自协方差
对于随机信号x（t），自相关函数为：
（9-13）
式中，
为时移。
若去掉x（t）的均值部分，则相应的自相关函数称为自协方差，即：
（9-14）
（9-15）
对于离散随机信号序列，x（n）的自相关函数和自协方差为：
（9-16）
（9-17）
式中，m为延迟。
9.2.2互相关函数和互协方差
对于两个不同随机信号x（t）,y（t），互相关函数为：
（9-18）
式中，
为时移。
互协方差为：
（9-19）
对于互协方差有：
（9-20）
对于离散随机序列x（n）和y（n），互相关函数和互协方差为：
（9-21）
（9-22）
注意到上面的公式中要求
或
，而在工程中信号是有限长的，因此只能得到相关函数和协方差的估计值，当t或N足够长时，估计值能精确地逼近真实值。
9.2.3计算相关函数和协方差的MATLAB函数
MATLAB信号处理工具箱提供了计算随机信号相关函数xcorr。
函数xcorr用于计算随机序列自相关和互相关函数。
调用格式为：
[c,lags]=xcorr（x,y[,maxlags,’option’]）
式中，x,y为两个独立的随机信号序列，长度均为N；c为x,y的互相关估计;lags为相关估计c的序号向量，其范围为[-maxlags:
maxlags]。
option缺省或’none’时，函数xcorr按下式执行非归一化的相关：
（9-23）
上角加星号的变量表示其共轭。
option为’biased’,计算有偏互相关函数估计
（9-24）
option为’unbiased’,计算无偏互相关函数估计
（9-25）
option为’coeff’,序列归一化，使零延迟的自相关函数为1。
maxlags为x和y的最大延迟,该项缺省时，函数返回值c的长度为2N-1;该项不缺省时，函数返回值c的长度为2*maxlags+1。
该函数也可用于求一个随机信号序列x（n）的自相关函数，调用格式为：
c=xcorr（x,maxlags）
MATLAB信号处理工具箱还提供了计算协方差函数xcov。
若x为一个向量，调用c=xcov（x）函数，则输出c为x序列的方差；若x为一个矩阵，则c=xcov（x）输出矩阵x的协方差矩阵，其对角线元素为每一列数据的方差，非对角线元素表示对应行和对应列的协方差。
调用c=xcov（x,y），其中x,y具有相同的长度，则c为式（9-22）所对应的互协方差。
【例9-2】求带有白噪声干扰的频率为10Hz的正弦信号和白噪声信号的自相关函数并进行比较。
%Samp9_2
clf;N=1000;Fs=500;%数据长度和采样频率
n=0:
N-1;t=n/Fs;%时间序列
Lag=100;%延迟样点数
randn（'state',0）;%设置产生随机数的初始状态
x=sin（2*pi*10*t）+0.6*randn（1,length（t））;%原始信号
[c,lags]=xcorr（x,Lag,'unbiased'）;%对原始信号进行无偏自相关估计
subplot（2,2,1）,plot（t,x）;%绘原始信号x
xlabel（'时间/s'）;ylabel（'x（t）'）;title（'带噪声周期信号'）;gridon;
subplot（2,2,2）;plot（lags/Fs,c）;%绘x信号自相关,lags/Fs为时间序列
xlabel（'时间/s'）;ylabel（'Rx（t）'）;
title（'带噪声周期信号的自相关'）;gridon;
%信号x1
x1=randn（1,length（x））;%产生一与x长度一致的随机信号
[c,lags]=xcorr（x1,Lag,'unbiased'）;%求随机信号x1的无偏自相关
subplot（2,2,3）,plot（t,x1）;%绘制随机信号x1
xlabel（'时间/s'）;ylabel（'x1（t）'）;title（'噪声信号'）;
gridon;
subplot（2,2,4）;plot（lags/Fs,c）;%绘制随机信号x1的无偏自相关
xlabel（'时间/s'）;ylabel（'Rx1（t）'）;
title（'噪声信号的自相关'）;
gridon
图9-1例9-2求得带有白噪声的周期信号和噪声信号的自相关
程序运行结果见图9-1。
可以看出，含有周期成分和噪声成分的自相关函数在
时具有最大值。
且在
较大时仍具有明显的周期性，其频率和周期信号的频率相同；而不含周期成分的纯噪声信号在
时具有最大值。
且在
稍大时很快衰减至零。
自相关的这一性质被用来识别随机信号中是否含有周期成分并用于确定所含周期成分的频率。
【例9-3】已知两个周期信号：
，其中f=10Hz，求互相关函数
。
%Samp9_3
clf;N=1000;Fs=500;%数据长度和采样频率
n=0:
N-1;t=n/Fs;%时间序列
Lag=200;%最大延迟单位
x=sin（2*pi*10*t）;%周期信号x
y=0.5*sin（2*pi*10*t+90*pi/180）;%与x有90o相移的信号y
[c,lags]=xcorr（x,y,Lag,'unbiased'）;%求无偏互相关
subplot（2,1,1）;plot（t,x,'r'）;%绘制x信号
holdon;plot（t,y,':
'）;%在同一幅图中绘y信号
legend（'x信号','y信号'）
xlabel（'时间/s'）;ylabel（'x（t）y（t）'）;
title（'原始信号'）;gridon;
holdoff
subplot（2,1,2）,plot（lags/Fs,c,'r'）;%绘制x,y的互相关
xlabel（'时间/s'）;ylabel（'Rxy（t）'）;
title（'信号x和y的相关'）;gridon
图9-2例9-3所给的原始信号及其互相关
程序运行结果为图9-2。
由结果可见，
也是周期为10Hz的余弦信号（注意两个图形其中的时间尺度有改变），其幅值为0.25，初始相位为90o。
由此可得到互相关函数的一个重要特性：
两个均值为零且具有相同频率的周期信号，其互相关函数
保留原信号频率、相位差和幅值信息。
下面例子给出互相关函数的另一个工程应用实例。
若信号x（t）和信号y（t）是同类信号且有时移，用互相关函数可以准确地计算出两个信号时移大小。
这种特性使得互相关函数广泛地应用于工程测量技术中。
【例9-4】两个sinc信号有0.2s的时移，用互相关函数计算时移的大小。
%Samp9_4
clf
N=1000;n=0:
N-1;Fs=500;t=n/Fs;%数据个数采样频率和时间序列
Lag=200;%最大延迟单位数
x1=90*sinc（pi*（n-0.1*Fs））;%第一个原始信号，延迟0.1s
y1=50*sinc（pi*（n-0.3*Fs））;%第二个原始信号,延迟0.3s
[c,lags]=xcorr（x1,y1,Lag,'unbiased'）;%计算两个函数的互相关
subplot（2,1,1）,plot（t,x1,'r'）;%绘第一个信号
holdon;plot（t,y1,'b:
'）;%在同一幅图中绘第二个信号
legend（'信号x','信号y'）;%绘制图例
xlabel（'时间/s'）;ylabel（'x（t）y（t）'）;
title（'信号x和y'）;holdoff
subplot（2,1,2）,plot（lags/Fs,c,'r'）;%绘制互相关信号
xlabel（'时间/s'）;ylabel（'Rxy（t）'）;
title（'信号x和y的相关'）;
程序运行结果见图9-3。
可以清楚地看到第二个信号相对于第一个信号延迟了0.2s，即在-0.2s处出现相关极大值。
因此可以采用该项技术检测延迟信号。
图9-3例9-4两个sinc信号的互相关
9.3功率谱估计
功率谱估计的目的是根据有限数据给出信号、随机过程的频率成分分布的描述。
上节所述的相关分析是时域内在噪声背景下提取有用信息的途径，而功率谱是频域内提取淹没在噪声中有用信息的分析方法。
9.3.1功率谱密度
假如随机信号x（t）的自相关函数为
，
的Fourier变换为：
（9-26）
则定义Sx（f）为x（t）的自功率谱密度或称为自功率谱。
因为Sx（f）可解释为x（t）的平均功率相对于频率的分布函数。
自功率谱Sx（f）包含
的全部信息。
如果噪声信号中含有某种频率成分，可以从自功率谱中看出。
若随机信号x（t）的自功率谱为Sx（f），则根据Fourier逆变换可得：
（9-27）
和自功率谱类似，两个随机信号x（t）和y（t）互相关的频率特性可用互功率谱密度来描述。
互功率谱密度和互相关函数也是一Fourier变换对。
（9-28）
（9-29）
对于离散随机序列x（n），自功率谱密度Sx（f）和自相关函数Rx（m）的关系为：
（9-30）
这里，Ts为数据采样间隔。
对于离散随机序列x（n）和y（n），互功率谱密度Sxy（f）和互相关函数Rxy（m）关系为：
（9-31）
且有
（9-32）
实际工程中随机序列长度均为有限长，因此利用有限长随机序列计算的自功率谱密度和互功率谱密度只是真实值的一种估计。
功率谱估计方法一般可分为参数估计和非参数估计两类。
MATLAB信号处理工具箱介绍的功率谱密度非参数估计方法有WELCH法、MTM法和MUSIC法；参数估计方法有MEM法，下面先介绍求取功率谱估计的基本方法，然后逐一介绍信号处理工具箱中给出的上述方法。
9.3.2周期图法
1.基本方法
周期图法是直接将信号的采样数据x（n）进行Fourier变换求取功率谱密度估计的方法。
假定有限长随机信号序列为x（n）。
它的Fourier变换和功率谱密度估计
存在下面的关系：
（9-33）
式中，N为随机信号序列x（n）的长度。
在离散的频率点
有：
（9-34）
其中，FFT[x（n）]为对序列x（n）的Fourier变换，由于FFT[x（n）]的周期为N，求得的功率谱估计以N为周期，因此这种方法称为周期图法。
下面用例子说明如何采用这种方法进行功率谱估计。
【例9-5】用Fourier变换算法求信号
的功率谱。
其中f1=50Hz,f2=120Hz,w（t）为白噪声，采样频率Fs=1000Hz。
（1）信号长度N=256；
（2）信号长度N=1024。
%Samp9_5
clf;Fs=1000;
%第一种情况:
N=256
N=256;Nfft=256;%数据长度和FFT所用的数据长度
n=0:
N-1;t=n/Fs;%采用的时间序列
xn=sin（2*pi*50*t）+2*sin（2*pi*120*t）+randn（1,N）;%带有噪声的信号
Pxx=10*log10（abs（fft（xn,Nfft）.^2）/N）;%Fourier振幅谱平方的平均值，并转换为dB
f=（0:
length（Pxx）-1）*Fs/length（Pxx）;%给出频率序列
subplot（2,1,1）,plot（f,Pxx）;%绘制功率谱曲线
xlabel（'频率/Hz'）;ylabel（'功率谱/dB'）;
title（'周期图N=256'）;gridon;
%第二种情况:
N=1024
Fs=1000;N=1024;Nfft=1024;%采样频率、数据长度和FFT所用数据长度
n=0:
N-1;t=n/Fs;%时间序列
xn=sin（2*pi*50*t）+2*sin（2*pi*120*t）+randn（1,N）;%带有噪声原始信号
Pxx=10*log10（abs（fft（xn,Nfft）.^2）/N）;%Fourier振幅谱平方的平均值，并转换为dB
f=（0:
length（Pxx）-1）*Fs/length（Pxx）;%频率序列
subplot（2,1,2）,plot（f,Pxx）;%绘制功率谱曲线
xlabel（'频率/Hz'）;ylabel（'功率谱/dB'）;
title（'周期图N=1024'）;gridon
程序的运行结果为图9-4。
可以看出，在频率50Hz和120Hz处，功率谱有两个峰值，说明信号中含有50Hz和120Hz两种频率成分。
但功率谱密度在很大范围内波动，而且并没有因信号取样点数N的增加而有明显改进。
注意本例中500Hz为Nyquist频率，由第3章介绍的Fourier变换的知识可知，只观看500Hz之前的功率谱可以推知以后的功率谱。
在本章后面的例题中我们将只研究Nyquist频率之前的功率谱。
图9-4例9-5中含有噪声信号的功率谱
用有限长样本序列的Fourier变换来表示随机序列的功率谱，只是一种估计或近似，不可避免存在误差。
为了减少误差，使功率谱估计更加平滑，可采用分段平均周期图法（Bartlett法）、加窗平均周期图法（Welch法）等方法加以改进。
2.分段平均周期图法（Bartlett法）
将信号序列x（n），n=0,1,…,N-1，分成互不重叠的P个小段，每小段由m个采样值，则P*m=N。
对每个小段信号序列进行功率谱估计，然后再取平均作为整个序列x（n）的功率谱估计。
平均周期图法还可以对信号x（n）进行重叠分段，如按2:
1重叠分段，即前一段信号和后一段信号有一半是重叠的。
对每一小段信号序列进行功率谱估计，然后再取平均值作为整个序列x（n）的功率谱估计。
这两种方法都称为平均周期图法，一般后者比前者好。
【例9-6】对例9-5中的信号序列，用两种平均周期图法求信号的功率谱密度估计。
%Samp9_6
clf;Fs=1000;%数据采样点数
%运用信号不重叠分段估计功率谱
N=1024;Nsec=256;n=0:
N-1;t=n/Fs;%数据点数，分段间隔，时间序列
randn（'state',0）;%设置产生随机数的初始状态
xn=sin（2*pi*50*t）+2*sin（2*pi*120*t）+randn（1,N）;%带噪声的原始信号
pxx1=abs（fft（xn（1:
256）,Nsec）.^2）/Nsec;%第一段功率谱
pxx2=abs（fft（xn（257:
512）,Nsec）.^2）/Nsec;%第二段功率谱
pxx3=abs（fft（xn（515:
768）,Nsec）.^2）/Nsec;%第三段功率谱
pxx4=abs（fft（xn（769:
1024）,Nsec）.^2）/Nsec;%第四段功率谱
Pxx=10*log10（（pxx1+pxx2+pxx3+pxx4）/4）;%平均得到整个序列功率谱
f=（0:
length（Pxx）-1）*Fs/length（Pxx）;%给出功率谱对应的频率
subplot（2,1,1）,plot（f（1:
Nsec/2）,Pxx（1:
Nsec/2））;%绘制功率谱曲线
xlabel（'频率/Hz'）;ylabel（'功率谱/dB'）;
title（'平均周期图（无重叠）N=4*256'）;
gridon
%运用信号重叠分段估计功率谱
pxx1=abs（fft（xn（1:
256）,Nsec）.^2）/Nsec;%第一段功率谱
pxx2=abs（fft（xn（129:
384）,Nsec）.^2）/Nsec;%第二段功率谱
pxx3=abs（fft（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（fft（xn（769:
1024）,Nsec）.^2）/Nsec;%第七段功率谱
Pxx=10*log10（（pxx1+pxx2+pxx3+pxx4+pxx5+pxx6+pxx7）/7）;%功率谱平均并转化为dB
f=（0:
length（Pxx）-1）*Fs/length（Pxx）;%频率序列
subplot（2,1,2）,plot（f（1:
Nsec/2）,Pxx（1:
Nsec/2））;%绘制功率谱曲线
xlabel（'频率/Hz'）;ylabel（'功率谱/dB'）;
title（'平均周期图（重叠一半）N=1024'）;
gridon
图9-5例9-6中重叠和不重叠分段计算的功率谱比较
程序运行结果为图9-5，上图采用不重叠分段法的功率谱估计，下图为2:
1重叠分段的功率谱估计，可见后者估计曲线较为平滑。
与上例比较，平均周期图法功率谱估计具有明显效果（涨落曲线靠近0dB）。
3.加窗平均周期图法
加窗平均周期图法是对分段平均周期图法的改进。
在信号序列x（n）分段后，用非矩形窗口对每一小段信号序列进行预处理，再采用前述分段平均周期图法进行整个信号序列x（n）的功率谱估计。
由窗函数的基本知识（第7章）可知，采用合适的非矩形窗口对信号进行处理可减小“频谱泄露”，同时可增加频峰的宽度，从而提高频谱分辨率。
【例9-7】对例9-5中的信号序列，用加窗平均周期图法进行功率谱密度估计。
%Samp9_7
clf;Fs=1000;%采样频率
N=1024;Nsec=256;n=0:
N-1;t=n/Fs;%数据长度，分段数据长度、时间序列
w=hanning（256）';%采用的窗口数据
randn（'state',0）;%设置产生随机数的初始状态
xn=sin（2*pi*50*t）+2*sin（2*pi*120*t）+randn（1,N）;%带噪声的信号
%采用不重叠加窗方法的功率谱估计
pxx1=abs（fft（w.*xn（1:
256）,Nsec）.^2）/norm（w）^2;%第一段加窗振幅谱平方
pxx2=abs（fft（w.*xn（25
升级会员 