《MATLAB辅助现代工程数字信号处理》课件第6章

第6章平稳随机信号处理与分析6.1随机信号及其处理

6.2平稳随机信号的时域描述

6.3平稳随机信号的频域描述

6.4线性系统对随机信号的响应

6.5平稳随机信号的模型

6.6小结

6.1随机信号及其处理

6.1.1随机信号处理的发展历程

随机信号处理的发展可分为两个阶段:经典随机信号处理阶段和现代随机信号处理阶段。

第一阶段为经典随机信号理论和技术的生长、发展和成熟期。

第二阶段为现代随机信号处理理论与技术起步和大发展的时期。6.1.2随机信号及其特征描述

随机信号(或序列)是一个随机过程，在它的每个时间点上的取值都是随机的，可用一个随机变量表示。一个随机过程是一个随机试验所产生的随机变量依时序组合得到的序列。随机信号X(t)是依赖时间t的随机变量，可以用描述随机变量的方法来描述随机信号。

当t在时间轴上取值t1，t2，…，tm时，可得到m个随机变量X(t1)，X(t2)，…，X(tm)，该随机信号可利用m维的概率分布函数(或概率密度)来描述:(6.1)

1.均值(一阶矩)

离散随机信号X(n)的所有样本函数，在同一时刻取值的统计平均值称为集平均，简称均值。即

(6.2)

2.方差(二阶矩)方差用于说明随机信号各可能值对其平均值的偏离程度，是随机信号在均值上下起伏变化的一种度量，它定义为可能值偏离其均值平方的数学期望。方差表示为(6.3)

方差的平方根称为均方差或标准差，即(6.4)

3.均方值

均方值描述了离散随机信号的强度或功率。即(6.5)

均方值与离散随机信号的均值和方差的关系为(6.6)4.自相关函数自相关函数可表示为(6.7)

5.自协方差函数

自协方差函数可表示为

(6.8)随机信号的自相关函数φX(n1，n2)描述了信号X(n)在n1、n2这两个时刻的相互关系，是一个重要的统计量。对于两个随机信号X(n)、Y(n)，其互相关函数和互协方差函数分别为

(1)互相关函数:(6.9)(2)互协方差函数:(6.10)

6.2平稳随机信号的时域描述

6.2.1平稳随机信号的数字特征

平稳随机信号是一类重要的随机信号。

一个离散随机信号X(n)，如果其均值与时间n无关，其自相关函数φX(n1，n2)与n1和n2的选取无关，而仅与n1和n2之差有关，那么称X(n)为广义平稳随机信号。狭义的平稳随机信号是指概率特性不随时间的平移而变化(与时间基准点无关)的随机信号。只有当X(n)是高斯随机过程时，才是狭义的平稳。

对于平稳随机信号，其均值、方差及均方值描述如下。

1.均值

均值可表示为(6.11)对于有限长平稳随机信号序列，为计算其均值估计，均值可表示为(6.12)

MATLAB提供了函数mean来计算随机离散信号的均值，调用格式为

y=mean(x)其中，x为离散随机序列，y为其均值。

2.方差

方差可表示为

(6.13)

MATLAB提供了函数std来计算随机离散信号的均方差(标准差)，调用格式为

y=std(x)

y=std(x，flag)其中，x为离散随机序列;y为方差;flag为控制符，用来控制计算均方差的算法，当flag=0时，计算无偏均方差(6.14)当flag=1时，计算有偏均方差为(6.15)

3.均方值

均方值可表示为(6.16)

【例6.1】

利用MATLAB编程计算长度N=100000的正态分布高斯随机噪声的均值、均方值、均方根值、方差以及均方差。

MATLAB程序如下:

%MATLABPROGRAM6-1

N=10^5;

x=randn(1，N);%产生正态分布高斯随机噪声

disp(′均值′);x_mean=mean(x)

disp(′均方值′);x_2p=x*x′/N

disp(′均方根值′);x_sq=sqrt(x_2p)

disp(′均方差′);x_std=std(x，1)

disp(′方差′);x_d=x_std.*x_std

程序运行结果为

均值

x_mean=0.0019

均方值

x_2p=0.9956

均方根值

x_sq=0.9978

均方差

x_std=0.9978

方差

x_d=0.9955

6.2.2相关函数和协方差

对于平稳随机信号来说，相关函数分为自相关函数和互相关函数;协方差分为自协方差和互协方差。

1.自相关函数

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

(6.17)2.自协方差函数对于随机信号x(t)，其自协方差函数为(6.18)3.互相关函数两个平稳随机信号X(n)、Y(n)的互相关函数为(6.19)

4.互协方差函数

两个平稳随机信号X(n)、Y(n)的互协方差函数为

(6.20)计算互相关函数的调用格式为

c=xcorr(x，y)

c=xcorr(x，y，′option′)

c=xcorr(x，y，maxlags，′option′)

[c，lags]=xcorr(x，y，maxlags，′option′)′option′为选择项:

(1)′biased′，计算有偏互相关函数估计:(6.21)

(2)′unbiased′，计算无偏互相关函数估计:(6.22)

(3)′coeff′，序列归一化，使零延迟的自相关函数为1。

(4)′none′，option缺省情况，此时，函数xcorr按下式执行非归一化计算行相关:(6.23)

【例6.2】

已知两个周期信号:x(t)=2sin(2πft)，y(t)=sin(2πft+85°)。其中，f＝8Hz，求两信号的互相关函数。

MATLAB程序如下:

%MATLABPROGRAM6-2

clf;

N=500;

n=[0:N－1];

fs=500;

t=n/fs;

lag=200;

x=2*sin(2*pi*8*t);

y=sin(2*pi*8*t+85*pi/180);

[c，lags]=xcorr(x，y，lag，′unbiased′);

subplot(221);

plot(t，x，′r－′，t，y，′b:′);

xlabel(′t′);ylabel(′x(t)y(t)′);title(′原始随机序列′);

legend(′x(t)′，′y(t)′);

grid;

subplot(222)

plot(lags/fs，c);

xlabel(′t′);ylabel(′Rxy(t)′);title(′互相关函数′);

grid;

程序运行结果如图6.1所示。图6.1同频率周期信号的互相关函数

【例6.3】

已知正弦信号为x(t)=sin(2πft)，f=8Hz，求一白噪声加该信号的自相关函数以及白噪声的自相关函数。

MATLAB程序如下:

%MATLABPROGRAM6-3

clf;

N=1000;

n=[0:N－1];

fs=500;

t=n/fs;

lag=100;

%信号加噪声

x=sin(2*pi*8*t)+0.85*randn(1，length(t));

[c，lags]=xcorr(x，lag，′unbiased′);

subplot(221);

plot(t，x);xlabel(′t′);ylabel(′x(t)′);

title(′Originalsignalx′);grid;

subplot(222)

plot(lags/fs，c);

xlabel(′t′);ylabel(′Rx(t)′);

title(′Autocorrelation′);grid;

%噪声

d=randn(1，length(x));

[c，lags]=xcorr(d，lag，′unbiased′);

subplot(223);

plot(t，d);xlabel(′t′);ylabel(′d(t)′);

title(′Noisesignald′);grid;

subplot(224)

plot(lags/fs，c);

xlabel(′t′);ylabel(′Rd(t)′);

title(′Autocorrelation′);grid;

程序运行结果如图6.2所示。图6.2不同信号的自相关函数

【例6.4】

已知正弦信号为x(t)=sin(2πft)，归一化频率f=0.1Hz，利用MATLAB编程产生两个正弦加白噪声序列，求取两序列的协方差函数。

MATLAB程序如下:

%MATLABPROGRAM6-4

clf;

N=256;t=[0:N-1];

f=0.1;a1=5;a2=3;

Mlag=N/4;

%产生两个正弦加白噪声的数据

x=a1*sin(2*pi*f*t)+2*randn(1，N);

y=a2*sin(2*pi*f*t)+randn(1，N);

subplot(311);

plot(x(1:N/2));grid;

subplot(312);

plot(y(1:N/2));grid;

%求这两个数据向量的协方差函数

covxy=xcov(x，y，Mlag，′biased′);

subplot(3，1，3);

plot((－Mlag:1:Mlag)，covxy);

程序运行结果如图6.3所示。图6.3两序列的互协方差函数6.2.3平稳随机信号的各态遍历性

对一平稳随机信号X(n)，如果它的所有样本函数在某一固定时刻的一阶和二阶统计特性与单一样本函数在长时间内的统计特性一致，则称X(n)为各态遍历信号。单一样本函数随时间变化的过程可以包括该信号所有样本函数的取值遍历。

设x(n)是各态遍历信号X(n)的一个样本函数，对X(n)的数字特征可定义如下:

(6.24)(6.25)

【例6.5】

讨论随机相位正弦序列X(n)=Asin(2πfnTs+

φ)的各态遍历性。

解对X(n)=Asin(2πfnTs+j)，其单一的时间样本x(n)=

Asin(2πfnTs+j)，j为一常数，对X(n)作时间平均，显然:

由于上式是n对求和，故求和号中的第一项与n无关，而第二项应等于零，所以

6.3平稳随机信号的频域描述

功率谱密度P(ω)用来描述离散随机信号的功率在频域上的分布情况，反映了单位频带信号功率的大小，是频率的函数。设φx(m)、φxy(m)分别为离散随机信号x(n)的自相关函数以及x(n)和y(n)的互相关函数，则有(6.26)(6.27)

Px(ω)称为离散随机信号的自功率谱;Pxy(ω)称为离散随机信号x(n)和y(n)的互功率谱。假设离散随机信号的功率是有限的，那么其功率谱密度的反变换必然存在，为其相关函数，即(6.28)而(6.29)

E[|x(n)|2]代表信号的平均功率，这说明Px(ω)在－π≤ω≤π频域内的积分面积正比于信号的平均功率。因此，Px(ω)是x(n)的平均功率密度，Px(ω)又称为功率谱密度。一个平稳随机序列u(n)，如果其功率谱Pu(ejω)在|ω|≤π的范围内始终为一常数，如σu2，则称该序列为白噪声序列，其自相关函数为(6.30)

【例6.6】一白噪声序列分别通过一个低通滤波器和一个高通滤波器。低通滤波器性能指标为f=[00.60.71]，a=[1100]，w1=[1，10];高通滤波器指标为f=[00.50.71]，a=[0011]，w2=[101];fs=1。两个滤波器的输出分别是x和y，求它们的自功率谱和互功率谱。

MATLAB程序如下:

%MATLABPROGRAM6-6

clc;

%构造低通滤波系数b1和高通滤波序列b2

Fs=1;N=1024;

a=[1，1，0，0];f=[0，0.6，0.7，1];

weigh=[1，10];

b1=remez(42，f，a，weigh);

a=[0，0，1，1];

f=[0，0.5，0.7，1];

weigh=[10，1];

b2=remez(42，f，a，weigh);

[h1，w]=freqz(b1，1，256，1);

h1=abs(h1);

h1=20*log10(h1);

subplot(331);plot(w，h1);grid;

[h2，w]=freqz(b2，1，256，1);

h2=abs(h2);

h2=20*log10(h2);

subplot(334);plot(w，h2);grid;

%将高斯白噪通过两个滤波器，并分别计算自功率谱

r=randn(16384，1);

x1=filter(b1，1，r);

x=x1(50:50+N);

[xpsd，F]=psd(x，N/4，Fs);

xpsd=10*log10(xpsd);

subplot(332);plot(F，xpsd);grid;

y1=filter(b2，1，r);

y=y1(50:50+N);

[ypsd，F]=psd(y，N/4，Fs);

ypsd=10*log10(ypsd);

subplot(333);plot(F，ypsd);grid;

%估计x和y的互功率谱

[pxy，w]=csd(x，y，N/4，Fs，hamming(N/4)，0，′mean′);

pxy=20*log10(abs(pxy));

subplot(335);

plot(w，pxy);grid;

程序运行结果如图6.4所示。图6.4离散随机序列的功率谱6.4线性系统对随机信号的响应

设x(n)为一平稳随机信号，它通过一线性时不变系统H(z)后，输出y(n):

(6.31)因此，y(n)也是平稳随机的。若x(n)是确定性信号，则Y(ejω)=X(ejω)H(ejω)

(6.32)由于随机信号不存在傅里叶变换，因此需从相关函数和功率谱的角度来研究随机信号通过线性系统的行为。现假设x(n)是实信号，则y(n)也是实的。y(n)的均值为(6.33)因为x(n)是平稳随机过程，有(6.34)故(6.35)即当mx是与时间无关的常数时，my也是与时间无关的常数。6.4.1自相关函数及自功率谱

假设输出y(n)是平稳的，则y(n)的自相关函数φy(m)为(6.36)自功率谱密度为(6.37)6.4.2互相关函数和互功率谱

下面探论关于线性非时变系统的输入和输出之间的互相关函数φxy(m)。由定义可知

(6.38)

该式又称为输入-输出互相关定理。输出自相关函数为(6.39)设mx=0(自相关函数的Z变换存在)，将式(6.38)与式(6.39)转换到z域有:(6.40)(6.41)用功率谱表示为(6.42)(6.43)当输入为白噪声时，其功率谱密度Px(ω)为常数，按式(6.26)有(6.44)式(6.44)表明这个常数是σx2，故在白噪声情况下有(6.45)(6.46)式(6.45)说明白噪声的功率在频率轴上的分布密度处处相同(等于σx2)，并且它就等于输入信号的平均功率。将式(6.46)代入式(6.38)，得(6.47)代入式(6.42)，得(6.48)式(6.47)与式(6.48)说明由白噪声激励的线性非时变系统，其输入、输出互相关函数正比于系统的冲激响应h(m)，而其输入、输出的互功率谱正比于系统的频率响应H(ejω)。

MATLAB信号处理工具箱提供了函数tfe来估计LSI系统的频率响应，调用格式为

[H，F]=tfe(x,y,Nfft,Fs,windows,noverlap,dflag)

【例6.7】x是一均匀分布的白噪声，它通过一个低通滤波器的输出为y。该滤波器的性能指标为:频率f=[00.50.71]，幅值a=[1100]，权重w=[110];归一化采样频率fs=1Hz。试估计该系统的频率响应。

MATLAB程序如下:

%MATLABPROGRAM6-7

clc;

%求已知系统的幅频响应

Fs=1;N=256;

a=[1，1，0，0];f=[0，0.6，0.7，1];

weigh=[1，10];

b=remez(42，f，a，weigh);

[h，w]=freqz(b，1，N，1);

h=20*log10(abs(h));

subplot(221);plot(w，h);grid;

r=rand(4096，1);

x=filter(b，1，r);

M=N*2;

%由函数tfe辨识出的系统的幅频响应

[H，w]=tfe(r，x，M，Fs，hamming(M)，0，′mean′);

H=20*log10(abs(H));

subplot(222);

plot(w，H);grid;

程序运行结果如图6.5所示。由图可知，函数tfe辨识出的幅频响应的通带有“毛刺”，但大体的波形还是一致的。图6.5函数tfe实现系统辨识

【例6.8】

两个sinc信号:x(t)=80sinc(2πf(t－0.1));y(t)=40sinc(2πf(t－0.3))，其中，f=250Hz，N=1000，fs=

500Hz。试用互相关函数计算两信号时移的大小。

MATLAB程序如下:

%MATLABPROGRAM6-8

clf

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

f=250;

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

x=90*sinc(pi*(n-0.1*Fs));

%第一个原始信号，延迟0.1s

y=50*sinc(pi*(n-0.3*Fs));

%第二个原始信号，延迟0.3s

[c，lags]=xcorr(x,y,Lag,'unbiased');%计算两个函数的互相关

subplot(2，1，1)，plot(t，x，′r′);%绘第一个信号

holdon;plot(t，y，′b:′);%在同一幅图中绘第二个信号

legend(′x(t)′，′y(t)′);%绘制图例

xlabel(′t/s′);ylabel(′x(t)y(t)′);

title(′原始信号′);holdoff

subplot(2，1，2)，plot(lags/Fs，c，′r′);%绘制互相关信号

xlabel(′t/s′);ylabel(′Rxy(t)′);

title(′信号x和y的相关′);

程序运行结果如图6.6所示。由图可知，第二个信号y(t)相对于第一个信号x(t)延