《MATLAB輔助現(xiàn)代工程數(shù)字信號處理》課件第6章

第6章平穩(wěn)隨機信號處理與分析6.1隨機信號及其處理

6.2平穩(wěn)隨機信號的時域描述

6.3平穩(wěn)隨機信號的頻域描述

6.4線性系統(tǒng)對隨機信號的響應

6.5平穩(wěn)隨機信號的模型

6.6小結

6.1隨機信號及其處理

6.1.1隨機信號處理的發(fā)展歷程

隨機信號處理的發(fā)展可分為兩個階段:經(jīng)典隨機信號處理階段和現(xiàn)代隨機信號處理階段。

第一階段為經(jīng)典隨機信號理論和技術的生長、發(fā)展和成熟期。

第二階段為現(xiàn)代隨機信號處理理論與技術起步和大發(fā)展的時期。6.1.2隨機信號及其特征描述

隨機信號(或序列)是一個隨機過程，在它的每個時間點上的取值都是隨機的，可用一個隨機變量表示。一個隨機過程是一個隨機試驗所產(chǎn)生的隨機變量依時序組合得到的序列。隨機信號X(t)是依賴時間t的隨機變量，可以用描述隨機變量的方法來描述隨機信號。

當t在時間軸上取值t1，t2，…，tm時，可得到m個隨機變量X(t1)，X(t2)，…，X(tm)，該隨機信號可利用m維的概率分布函數(shù)(或概率密度)來描述:(6.1)

1.均值(一階矩)

離散隨機信號X(n)的所有樣本函數(shù)，在同一時刻取值的統(tǒng)計平均值稱為集平均，簡稱均值。即

(6.2)

2.方差(二階矩)方差用于說明隨機信號各可能值對其平均值的偏離程度，是隨機信號在均值上下起伏變化的一種度量，它定義為可能值偏離其均值平方的數(shù)學期望。方差表示為(6.3)

方差的平方根稱為均方差或標準差，即(6.4)

3.均方值

均方值描述了離散隨機信號的強度或功率。即(6.5)

均方值與離散隨機信號的均值和方差的關系為(6.6)4.自相關函數(shù)自相關函數(shù)可表示為(6.7)

5.自協(xié)方差函數(shù)

自協(xié)方差函數(shù)可表示為

(6.8)隨機信號的自相關函數(shù)φX(n1，n2)描述了信號X(n)在n1、n2這兩個時刻的相互關系，是一個重要的統(tǒng)計量。對于兩個隨機信號X(n)、Y(n)，其互相關函數(shù)和互協(xié)方差函數(shù)分別為

(1)互相關函數(shù):(6.9)(2)互協(xié)方差函數(shù):(6.10)

6.2平穩(wěn)隨機信號的時域描述

6.2.1平穩(wěn)隨機信號的數(shù)字特征

平穩(wěn)隨機信號是一類重要的隨機信號。

一個離散隨機信號X(n)，如果其均值與時間n無關，其自相關函數(shù)φX(n1，n2)與n1和n2的選取無關，而僅與n1和n2之差有關，那么稱X(n)為廣義平穩(wěn)隨機信號。狹義的平穩(wěn)隨機信號是指概率特性不隨時間的平移而變化(與時間基準點無關)的隨機信號。只有當X(n)是高斯隨機過程時，才是狹義的平穩(wěn)。

對于平穩(wěn)隨機信號，其均值、方差及均方值描述如下。

1.均值

均值可表示為(6.11)對于有限長平穩(wěn)隨機信號序列，為計算其均值估計，均值可表示為(6.12)

MATLAB提供了函數(shù)mean來計算隨機離散信號的均值，調(diào)用格式為

y=mean(x)其中，x為離散隨機序列，y為其均值。

2.方差

方差可表示為

(6.13)

MATLAB提供了函數(shù)std來計算隨機離散信號的均方差(標準差)，調(diào)用格式為

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的正態(tài)分布高斯隨機噪聲的均值、均方值、均方根值、方差以及均方差。

MATLAB程序如下:

%MATLABPROGRAM6-1

N=10^5;

x=randn(1，N);%產(chǎn)生正態(tài)分布高斯隨機噪聲

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相關函數(shù)和協(xié)方差

對于平穩(wěn)隨機信號來說，相關函數(shù)分為自相關函數(shù)和互相關函數(shù);協(xié)方差分為自協(xié)方差和互協(xié)方差。

1.自相關函數(shù)

對于隨機信號x(t)，自相關函數(shù)為

(6.17)2.自協(xié)方差函數(shù)對于隨機信號x(t)，其自協(xié)方差函數(shù)為(6.18)3.互相關函數(shù)兩個平穩(wěn)隨機信號X(n)、Y(n)的互相關函數(shù)為(6.19)

4.互協(xié)方差函數(shù)

兩個平穩(wěn)隨機信號X(n)、Y(n)的互協(xié)方差函數(shù)為

(6.20)計算互相關函數(shù)的調(diào)用格式為

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′，計算有偏互相關函數(shù)估計:(6.21)

(2)′unbiased′，計算無偏互相關函數(shù)估計:(6.22)

(3)′coeff′，序列歸一化，使零延遲的自相關函數(shù)為1。

(4)′none′，option缺省情況，此時，函數(shù)xcorr按下式執(zhí)行非歸一化計算行相關:(6.23)

【例6.2】

已知兩個周期信號:x(t)=2sin(2πft)，y(t)=sin(2πft+85°)。其中，f＝8Hz，求兩信號的互相關函數(shù)。

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(′互相關函數(shù)′);

grid;

程序運行結果如圖6.1所示。圖6.1同頻率周期信號的互相關函數(shù)

【例6.3】

已知正弦信號為x(t)=sin(2πft)，f=8Hz，求一白噪聲加該信號的自相關函數(shù)以及白噪聲的自相關函數(shù)。

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不同信號的自相關函數(shù)

【例6.4】

已知正弦信號為x(t)=sin(2πft)，歸一化頻率f=0.1Hz，利用MATLAB編程產(chǎn)生兩個正弦加白噪聲序列，求取兩序列的協(xié)方差函數(shù)。

MATLAB程序如下:

%MATLABPROGRAM6-4

clf;

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

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

Mlag=N/4;

%產(chǎn)生兩個正弦加白噪聲的數(shù)據(jù)

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;

%求這兩個數(shù)據(jù)向量的協(xié)方差函數(shù)

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

subplot(3，1，3);

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

程序運行結果如圖6.3所示。圖6.3兩序列的互協(xié)方差函數(shù)6.2.3平穩(wěn)隨機信號的各態(tài)遍歷性

對一平穩(wěn)隨機信號X(n)，如果它的所有樣本函數(shù)在某一固定時刻的一階和二階統(tǒng)計特性與單一樣本函數(shù)在長時間內(nèi)的統(tǒng)計特性一致，則稱X(n)為各態(tài)遍歷信號。單一樣本函數(shù)隨時間變化的過程可以包括該信號所有樣本函數(shù)的取值遍歷。

設x(n)是各態(tài)遍歷信號X(n)的一個樣本函數(shù)，對X(n)的數(shù)字特征可定義如下:

(6.24)(6.25)

【例6.5】

討論隨機相位正弦序列X(n)=Asin(2πfnTs+

φ)的各態(tài)遍歷性。

解對X(n)=Asin(2πfnTs+j)，其單一的時間樣本x(n)=

Asin(2πfnTs+j)，j為一常數(shù)，對X(n)作時間平均，顯然:

由于上式是n對求和，故求和號中的第一項與n無關，而第二項應等于零，所以

6.3平穩(wěn)隨機信號的頻域描述

功率譜密度P(ω)用來描述離散隨機信號的功率在頻域上的分布情況，反映了單位頻帶信號功率的大小，是頻率的函數(shù)。設φx(m)、φxy(m)分別為離散隨機信號x(n)的自相關函數(shù)以及x(n)和y(n)的互相關函數(shù)，則有(6.26)(6.27)

Px(ω)稱為離散隨機信號的自功率譜;Pxy(ω)稱為離散隨機信號x(n)和y(n)的互功率譜。假設離散隨機信號的功率是有限的，那么其功率譜密度的反變換必然存在，為其相關函數(shù)，即(6.28)而(6.29)

E[|x(n)|2]代表信號的平均功率，這說明Px(ω)在－π≤ω≤π頻域內(nèi)的積分面積正比于信號的平均功率。因此，Px(ω)是x(n)的平均功率密度，Px(ω)又稱為功率譜密度。一個平穩(wěn)隨機序列u(n)，如果其功率譜Pu(ejω)在|ω|≤π的范圍內(nèi)始終為一常數(shù)，如σu2，則稱該序列為白噪聲序列，其自相關函數(shù)為(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;

%構造低通濾波系數(shù)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(xiàn)]=psd(x，N/4，F(xiàn)s);

xpsd=10*log10(xpsd);

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

y1=filter(b2，1，r);

y=y1(50:50+N);

[ypsd，F(xiàn)]=psd(y，N/4，F(xiàn)s);

ypsd=10*log10(ypsd);

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

%估計x和y的互功率譜

[pxy，w]=csd(x，y，N/4，F(xiàn)s，hamming(N/4)，0，′mean′);

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

subplot(335);

plot(w，pxy);grid;

程序運行結果如圖6.4所示。圖6.4離散隨機序列的功率譜6.4線性系統(tǒng)對隨機信號的響應

設x(n)為一平穩(wěn)隨機信號，它通過一線性時不變系統(tǒng)H(z)后，輸出y(n):

(6.31)因此，y(n)也是平穩(wěn)隨機的。若x(n)是確定性信號，則Y(ejω)=X(ejω)H(ejω)

(6.32)由于隨機信號不存在傅里葉變換，因此需從相關函數(shù)和功率譜的角度來研究隨機信號通過線性系統(tǒng)的行為?，F(xiàn)假設x(n)是實信號，則y(n)也是實的。y(n)的均值為(6.33)因為x(n)是平穩(wěn)隨機過程，有(6.34)故(6.35)即當mx是與時間無關的常數(shù)時，my也是與時間無關的常數(shù)。6.4.1自相關函數(shù)及自功率譜

假設輸出y(n)是平穩(wěn)的，則y(n)的自相關函數(shù)φy(m)為(6.36)自功率譜密度為(6.37)6.4.2互相關函數(shù)和互功率譜

下面探論關于線性非時變系統(tǒng)的輸入和輸出之間的互相關函數(shù)φxy(m)。由定義可知

(6.38)

該式又稱為輸入-輸出互相關定理。輸出自相關函數(shù)為(6.39)設mx=0(自相關函數(shù)的Z變換存在)，將式(6.38)與式(6.39)轉(zhuǎn)換到z域有:(6.40)(6.41)用功率譜表示為(6.42)(6.43)當輸入為白噪聲時，其功率譜密度Px(ω)為常數(shù)，按式(6.26)有(6.44)式(6.44)表明這個常數(shù)是σx2，故在白噪聲情況下有(6.45)(6.46)式(6.45)說明白噪聲的功率在頻率軸上的分布密度處處相同(等于σx2)，并且它就等于輸入信號的平均功率。將式(6.46)代入式(6.38)，得(6.47)代入式(6.42)，得(6.48)式(6.47)與式(6.48)說明由白噪聲激勵的線性非時變系統(tǒng)，其輸入、輸出互相關函數(shù)正比于系統(tǒng)的沖激響應h(m)，而其輸入、輸出的互功率譜正比于系統(tǒng)的頻率響應H(ejω)。

MATLAB信號處理工具箱提供了函數(shù)tfe來估計LSI系統(tǒng)的頻率響應，調(diào)用格式為

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

【例6.7】x是一均勻分布的白噪聲，它通過一個低通濾波器的輸出為y。該濾波器的性能指標為:頻率f=[00.50.71]，幅值a=[1100]，權重w=[110];歸一化采樣頻率fs=1Hz。試估計該系統(tǒng)的頻率響應。

MATLAB程序如下:

%MATLABPROGRAM6-7

clc;

%求已知系統(tǒng)的幅頻響應

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;

%由函數(shù)tfe辨識出的系統(tǒng)的幅頻響應

[H，w]=tfe(r，x，M，F(xiàn)s，hamming(M)，0，′mean′);

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

subplot(222);

plot(w，H);grid;

程序運行結果如圖6.5所示。由圖可知，函數(shù)tfe辨識出的幅頻響應的通帶有“毛刺”，但大體的波形還是一致的。圖6.5函數(shù)tfe實現(xiàn)系統(tǒng)辨識

【例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。試用互相關函數(shù)計算兩信號時移的大小。

MATLAB程序如下:

%MATLABPROGRAM6-8

clf

N=1000;n=[0:N－1];Fs=500;t=n/Fs;%數(shù)據(jù)個數(shù)采樣頻率和時間序列

f=250;

Lag=200;%最大延遲單位數(shù)

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');%計算兩個函數(shù)的互相關

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)延