Matlab仿真窄帶隨機過程

隨機過程數(shù)學(xué)建模分析

任何通信系統(tǒng)都有發(fā)送機和接收機，為了提高系統(tǒng)的可靠性，即輸出信噪比，

通常在接收機的輸入端接有一個帶通濾波器，信道內(nèi)的噪聲構(gòu)成了一個隨機過

程，經(jīng)過該帶通濾波器之后，則變成了窄帶隨機過程，因此，討論窄帶隨機過程

的規(guī)律是重要的。

一、窄帶隨機過程。

一個實平穩(wěn)隨機過程X(t)，若它的功率譜密度具有下述性質(zhì)：

中心頻率為ω，帶寬為△ω=2ω，當(dāng)△ω<<ω 時，就可認為滿足窄帶條件。

c0c

若隨機過程的功率譜滿足該條件則稱為窄帶隨機過程。若帶通濾波器的傳輸函數(shù)

滿足該條件則稱為窄帶濾波器。隨機過程通過窄帶濾波器傳輸之后變成窄帶隨機

過程。

圖1 為典型窄帶隨機過程的功率譜密度圖。若用一示波器來觀測次波形，則

可看到，它接近于一個正弦波，但此正弦波的幅度和相位都在緩慢地隨機變化，

圖2所示為窄帶隨機過程的一個樣本函數(shù)。

圖1 典型窄帶隨機過程的功率譜密度圖

圖2 窄帶隨機過程的一個樣本函數(shù)

二、窄帶隨機過程的數(shù)學(xué)表示

1、用包絡(luò)和相位的變化表示

由窄帶條件可知,窄帶過程是功率譜限制在ω附近的很窄范圍內(nèi)的一個隨機

c

過程,從示波器觀察(或由理論上可以推知):這個過程中的一個樣本函數(shù)(一個實現(xiàn))

的波形是一個頻率為?且幅度和相位都做緩慢變化的余弦波。

c

寫成包絡(luò)函數(shù)和隨機相位函數(shù)的形式：

X(t)=A(t)*cos[ωt+ Φ(t)]

c

其中:A(t)稱作X(t)的包絡(luò)函數(shù); Φ(t)稱作X(t)的隨機相位函數(shù)。包絡(luò)隨時間做

緩慢變化，看起來比較直觀，相位的變化，則看不出來。

2、萊斯（Rice）表示式

任何一個實平穩(wěn)隨機過程X(t)都可以表示為：

X(t)=A(t) cosωt-A(t) sinωt

ccSc

其中同相分量：

At＋ sinωt=LP[X(t) *2cosωt]

cccc

(t)= X(t) cosφt= X(t) cosω

正交分量：

A(t) = X(t)sinφt= cosωt— X(t) sinωt= LP[-X(t) *2sinωt]

Sccc

（LP[A]表示取A的低頻部分）。A(t)和A(t)都

cS

是實隨機過程，均值為0，方差等于X(t)的方差。

三、窄帶隨機過程仿真建模要求

1、用Matlab 編程仿真窄帶隨機信號：X(t)=(1+ A(t))*cos(ωt+φ)+n(t)。 其中

c

包絡(luò)A(t)頻率為1KHz，幅值為l V。載波頻率為：4KHz，幅值為l V，φ是一個

固定相位，n(t)為高斯白噪聲，采樣頻率設(shè)為16KHz。實際上，這是一個帶有載

波的雙邊帶調(diào)制信號。

2、計算窄帶隨機信號的均值、均方值、方差、概率密度、頻譜及功率譜密度、

相關(guān)函數(shù)，用圖示法來表示。

3、窄帶系統(tǒng)檢測框圖如圖3所示。

圖3 窄帶系統(tǒng)檢測框圖

4、低通濾波器設(shè)計：

低通濾波器技術(shù)要求：通帶截止頻率1KHz，阻帶截止頻率2KHz。過渡帶：

1KHz，阻帶衰減：>35DB，通帶衰減：<1DB，采樣頻率：≤44.1KHz

5、計算a點、b點、A(t)、A(t)、y(t)的均值、均方值、方差、頻譜及功率譜

cS

密度、相關(guān)函數(shù)，用圖示法來表示。

四、建模仿真過程及結(jié)果（程序見附件）

1、根據(jù)要求得到X(t)的表達式：

x= (l+a) .*cos (2*pi*4000*t+2) +noisy/10;

其中：noisy為高斯白噪聲，由wgn函數(shù)生成，

a=cos (2*pi*l000*t)，

均值：Ex=mean (x)，

方差：Dx=var (x)，

計算可得：X(t)的均值為0.0019，

X(t)的方差為0.7590。

如圖4所示，其中藍色線為X(t)一個樣本的時域波形，紅色點連成的線為X(t)

的均值，綠色點連成的線為X(t)的方差。

圖4 窄帶隨機信號時域波形

2、求X(t)的概率密度，方法是將最大最小區(qū)間分成14等份，然后分別計算各個

區(qū)間的個數(shù)，如圖02中柱形條所示，利用曲線擬合， 得到合適的概率密度函數(shù)。

為了得到光滑的曲線，利用了多項式擬合，經(jīng)過測試，9次擬合曲線比較符合要

求，獲得的曲線如圖5中曲線所示：

圖5 X(t)的概率分布密度函數(shù)

3、對X(t)進行頻譜分析，在Matlab中，利用fft函數(shù)可以很方便得求得X(t)的頻

譜，然后用abs和angle函數(shù)求得幅值和相位，畫出圖像如圖6所示：

圖6 X(t)的頻譜圖

4、求X(t)的自相關(guān)函數(shù)，用xcorr函數(shù)求出自相關(guān)序列，得到X(t)自相關(guān)函數(shù)的

時域波形，如圖7所示。

圖7 X(t)自相關(guān)函數(shù)的時域波形

5、對X(t)自相關(guān)函數(shù)進行fft變換，得到X(t)的功率譜密度，如圖8所示。

圖8 X(t)的功率譜密度

6、建立濾波器，建立一個巴特沃思濾波器，對產(chǎn)生的x(t)進行檢測。濾波器的幅

度譜和相位譜所示：

圖9 地通濾波器的幅度譜和相位譜

7、求A(t)的統(tǒng)計特性，A(t)為X(t) *2cosωt通過低通濾波器的信號，

ccc

A(t)的均值Eh = -0.4075 4（帶有直流分量），

c

A(t)的均方值是E2h =0.2458

c

A(t)的方差Dh = 0.0798

c

A(t)的波形如圖10、圖11所示：

c

圖10 A(t)的時域波形圖和頻譜圖

c

圖11 A(t)的自相關(guān)函數(shù)的時域波形圖和Ac(t)的功率譜密度

c

8、求A(t)的統(tǒng)計特性，A(t)為X(t) *2cosωt通過低通濾波器的信號，

SSc

A(t)的均值Eh =0.8972（帶有直流分量），

S

A(t)的均方值是E2h = 1.1565

S

A(t)的方差Dh = 0.3518

S

A(t)的波形如圖13、圖14所示：

S

圖13 A(t)的時域波形圖和頻譜圖

S

圖14 A(t)的自相關(guān)函數(shù)的時域波形圖和A(t)的功率譜密度

SS

9、求出Y(t)的統(tǒng)計特性，Y (t)=A(t) cosωt-At，

ccSc

(t) sinω

其統(tǒng)計特性如下

輸出信號Y(t)的均值Eh = -4.4011e-004s

輸出信號Y(t)的均方值E2h = 3.0280

輸出信號Y(t)的方差Dh = 3.0303

Y(t)的仿真圖形如圖15、圖16所示。

圖15 Y(t)的時域波形圖和頻譜圖

圖16 Y(t)的自相關(guān)函數(shù)的時域波形圖和Y(t)的功率譜密度

附件：

clc

fs=16000; %設(shè)定采樣頻率

N=1300;

n=0:N-1; %取的樣本點數(shù)

t=n/fs; %獲得以1/16000為時間間隔采樣序列

noisy=wgn(1,N,0); %產(chǎn)生高斯白噪聲

a=cos(2*pi*1000*t); %獲取A(t)的采樣點

x=(1+a).*cos(2*pi*4000*t+2)+noisy/10; %獲取x(t)的采樣點

%以t為橫坐標(biāo)畫出x(t)的時域圖型

figure(1); subplot(2,1,1); plot(n,x);

axis([0 140 -3 3]);xlabel('采樣點');ylabel('X(t)/V');title('窄帶隨機信號波形');grid on;

%求X(t)的統(tǒng)計特性 并畫出來

disp('X(t)的均值為'); Ex=mean(x); disp(Ex);%求X(t)均值

hold on; plot(n,Ex,'r.');

disp('X(t)的方差為');Dx=var(x); disp(Dx);%求x(t)方差

hold on; plot(n,Dx,'g.');

%畫出X(t)的概率分布函數(shù)

each=linspace(min(x),max(x),14); %將最大最小區(qū)間分成14等份,然后分別計算各個區(qū)間的個數(shù)

nr=hist(x,each); %計算各個區(qū)間的個數(shù)

nr=nr/length(x); %計算各個區(qū)間的個數(shù)歸一化

subplot(2,1,2); p=polyfit(each,nr,9); %畫出概率分布直方圖

bar(each,nr); %多項式擬合

hold on; plot(each,nr,'g')

eachi=-2:0.1:2;

nri=polyval(p,eachi);

plot(eachi,nri,'r')

axis tight;title('X(t)概率密度分布');xlabel('X(t)');ylabel('P(x)');grid on;

%對X(t)進行頻譜分析

Fx=fft(x,N); %對x(t)進行fft變換，在0~16000區(qū)間內(nèi)得到2N-1個頻率值

magn=abs(Fx); %求x(t)幅值

xangle=angle(Fx); %求X(t)相位

labelang=(0:length(x)-1)*16000/length(x); %在0~16000區(qū)間內(nèi)求橫坐標(biāo)刻度

figure(2); plot(labelang,magn*10); %在0~16000區(qū)間內(nèi)做頻譜和相位圖

axis([0 16000 -0.5 600]); xlabel('頻率/Hz');ylabel('幅值');title('X(t)頻譜圖');grid on;

%求X(t)的自相關(guān)函數(shù)

[c,lags]=xcorr(x,'coeff'); %求出自相關(guān)序列

figure(3); subplot(2,1,1); plot(lags/fs,c); %在時域內(nèi)畫自相關(guān)函數(shù)

axis tight; xlabel('T');ylabel('Rx(T)');title('X(t)的自相關(guān)函數(shù)');grid on;

%求X(t)的功率譜密度

long=length(c);

Sx=fft(c,long);

labelx=(0:long-1)*2*pi;

plot_magn=10*log10(abs(Sx));

subplot(2,1,2); plot(labelx,plot_magn); %畫功率譜密度

axis tight;xlabel('w');ylabel('Sx(w)');title('X(t)的功率譜密度');grid on;

%窄帶系統(tǒng)檢測

z1=2.*cos(2*pi*4000*t);

z2=-2.*sin(2*pi*4000*t);

Ac=z1.*x; %濾波后生成Ac(t)

As=z2.*x; %濾波后生成As(t)

y=Ac.*cos(2*pi*4000*t)-As.*sin(2*pi*4000*t);

%濾波器設(shè)計

f_p=1000;f_s=1600;R_p=1;R_s=35; %設(shè)定濾波器參數(shù); 通、阻帶截止頻率,通、阻帶衰減

Ws=2*f_s/fs;Wp=2*f_p/fs; %頻率歸一化

[n,Wn]=buttord(Wp,Ws,R_p,R_s); %采用巴特沃思濾波器

[b,a]=butter(n,Wn); %求得濾波器傳輸函數(shù)的多項式系數(shù)

figure(4);

[H,W]=freqz(b,a); %求得濾波器傳輸函數(shù)的幅頻特性

subplot(2,1,1); plot(W*fs/(2*pi),abs(H)); %在0～2pi區(qū)間內(nèi)作幅度譜

title('低通濾波器幅度譜'); grid on;

subplot(2,1,2); plot(W*fs/(2*pi),angle(H)); %在0～2pi區(qū)間內(nèi)作相位譜

title('低通濾波器相位譜'); grid on;

%求Ac(t)濾波后的統(tǒng)計特性

mc=filter(b,a,Ac); %上支路通過濾波器 Ac(t)

disp('Ac(t)的均值');Eh=mean(mc) %求Ac(t)的均值

disp('Ac(t)的均方值是');E2h=mc*mc'/N %求Ac(t)的均方值

disp('Ac(t)的方差');Dh=var(mc) %求Ac(t)的方差

%畫Ac(t)的時域波形

figure(6); subplot(2,1,1); n=0:N-1; plot(n,mc);

axis([0 300 -1 1]);xlabel('采樣點');ylabel('幅值');title('Ac(t)的時域波形');grid on;

%畫Ac(t)的頻譜圖

yc=fft(mc,length(mc)); %對Ac(t)進行fft變換

longc=length(yc); %求傅里葉變換后的序列長度

labelx=(0:longc-1)*16000/longc;

magnl=abs(yc); %求Ac(t)的幅值

subplot(2,1,2); plot(labelx,magnl); %畫Ac(t)的頻譜圖

axis tight; xlabel('頻率(Hz)'); ylabel('幅值'); title('Ac(t)頻譜圖'); grid on;

%求Ac(t)的自相關(guān)函數(shù)

[c1,lags1]=xcorr(mc,'coeff'); %求出Ac(t)的自相關(guān)序列

figure(7); subplot(2,1,1); plot(lags1/fs,c1); %在時域內(nèi)畫Ac(t)的自相關(guān)函數(shù)

xlabel('T');ylabel('Rx(T)');axis tight;

title('Ac(t)的自相關(guān)函數(shù)');

grid on;

%求Ac(t)的雙邊功率譜

Sac=fft(c1,length(c1)); %對Ac(t)的自相關(guān)函數(shù)進行傅里葉變換

magnc=abs(Sac); %求Ac(t)的雙邊功率譜幅值

long=length(Sac); %求傅里葉變換后的序列長度

labelc=(0:long-1)*16000/long;

subplot(2,1,2); plot(labelc,10*log10(magnc)); %畫Ac(t)的自相關(guān)函數(shù)頻譜 即為Ac(t)的雙邊功率譜

xlabel('頻率(Hz)');ylabel('功率譜(dbW)');axis tight;title('Ac(t)的雙邊功率譜');grid on;

%求得As(t)的統(tǒng)計特性

ms=filter(b,a,As); %對下支路信號進行濾波得As(t)

disp('As(t)的均值'); Eh=mean(ms) %求As(t)的均值

disp('As(t)的均方值是'); E2h=ms*ms'/N %求As(t)的均方值

disp('As(t)的方差'); Dh=var(ms) %求As(t)的方差

%作As(t)的時域波形

figure(8);subplot(2,1,1); n=0:N-1;plot(n,ms); %畫出As(t)的時域波形

axis([0 300 -0.5 2]); xlabel('采樣點');ylabel('幅值');title('As(t)的時域波形');grid on;

%對As(t)進行FFT變換并做頻譜圖

ys=fft(ms,length(ms)); %對As(t)進行fft變換

longs=length(ys); %求傅里葉變換后的序列長度

labelx=(0:longs-1)*16000/longs;

magn2=abs(ys); %求As(t)的幅值

subplot(2,1,2); plot(labelx,magn2); %畫出As(t)的頻譜圖

axis tight;xlabel('頻率(Hz)');ylabel('幅值');title('As(t)的頻譜圖');grid on;

%求As(t)的自相關(guān)函數(shù)

[c2,lags2]=xcorr(ms,'coeff'); %求出As(t)的自相關(guān)序列

figure(9);subplot(2,1,1);plot(lags2/fs,c2); %畫出As(t)自相關(guān)函數(shù)的時域波形

xlabel('T');ylabel('Rx(T)');axis tight;title('As(t)的的自相關(guān)函數(shù)');grid on;

%求As(t)的雙邊功率譜

Sas=fft(c2,length(c2)); %對As(t)的自相關(guān)函數(shù)進行傅里葉變換

magnc=abs(Sac); %求As(t)的雙邊功率譜幅值

long=length(Sas); %求傅里葉變換后的序列長度

labels=(0:long-1)*16000/long;

subplot(2,1,2); plot(labelc,10*log10(magnc)); %畫As(t)的自相關(guān)函數(shù)頻譜

xlabel('頻率(Hz)');ylabel('功率譜(dbW)');axis tight;title('As(t)的雙邊功率譜');

% 求y(t)的統(tǒng)計特性

disp('輸出信號Y(t)的均值');Eh=mean(y) %求輸出信號Y(t)的均值

disp('輸出信號Y(t)的均方值');E2h=y*y'/N %求輸出信號Y(t)的均方值

disp('輸出信號Y(t)的方差');Dh=var(y) %求輸出信號Y(t)的方差

%作輸出信號Y(t)的時域波形

figure(10); subplot(2,1,1);n=0:N-1;plot(n,y);

axis([0 150 -2 2]);xlabel('采樣點');ylabel('幅值');title('Y(t)的時域波形');grid on;

%進行FFT變換并做頻譜圖

yy=fft(y,length(y)); %對相加后的信號進行fft變換

longy=length(yy); %Y(t)傅里葉變換后的序列長度

labelx=(0:longy-1)*16000/longy;

magn3=abs(yy); %求Y(t)的幅值

subplot(2,1,2); plot(labelx,magn3); %做Y(t)的頻譜圖

axis tight;xlabel('頻率(Hz)');ylabel('幅值');title('Y(t)的頻譜圖');grid on;

%求輸出信號Y(t)的自相關(guān)函數(shù)

[c3,lags3]=xcorr(y,'coeff'); %求出Y(t)的自相關(guān)序列

figure(11); subplot(2,1,1); plot(lags3/fs,c3); %畫Y(t)自相關(guān)函數(shù)的時域波形

xlabel('T');ylabel('Rx(T)');axis tight;title('Y(t)的的自相關(guān)函數(shù)');grid on;

%求輸出信號Y(t)的雙邊功率譜

Sy=fft(c3,length(c3)); %對Y(t)的自相關(guān)函數(shù)進行傅里葉變換

magny=abs(Sy); %求Y(t)雙邊功率譜幅值

long=length(Sy);

labely=(0:long-1)*16000/long;

subplot(2,1,2); plot(labely,10*log10(magny)); %****畫Y(t)的功率譜密度

xlabel('頻率(Hz)');ylabel('功率譜(dbW)');axis tight;title('Y(t)的雙邊功率譜');grid on;