白噪聲通過(guò)LTI-的仿真

1、實(shí)驗(yàn)2 白噪聲通過(guò)LTI 的仿真1、實(shí)驗(yàn)?zāi)康牧私獍自肼曂ㄟ^(guò)LTI 系統(tǒng)的原理與處理方法，學(xué)會(huì)運(yùn)用Matlab 函數(shù)對(duì)隨機(jī)過(guò)程進(jìn)行均值、相關(guān)函數(shù)和功率譜的估計(jì),并且通過(guò)實(shí)驗(yàn)分析理論分析與實(shí)驗(yàn)結(jié)果之間的差異。2、實(shí)驗(yàn)原理假定一具有單位方差的抽樣序列X(n)的白噪聲隨機(jī)過(guò)程X(t)通過(guò)一脈沖響應(yīng)為 的線性濾波器，繪出輸入輸出信號(hào)的均值、方差、相關(guān)函數(shù)及功率譜密度。設(shè)系統(tǒng)沖激響應(yīng)為h(n)，傳遞函數(shù) ，或者用Z 變換，結(jié)果為 。輸入為 X (n)，輸出為 ，均值關(guān)系： ，若平穩(wěn)有， 自相關(guān)函數(shù)關(guān)系， ，當(dāng)是平穩(wěn)時(shí)候，有 題目中假設(shè)為白噪聲，可以根據(jù)白噪聲的性質(zhì)進(jìn)行理論計(jì)算。白噪聲的自相關(guān)函數(shù)， 這里，

2、假設(shè)的是零均值和單位方差，于是，而 對(duì)應(yīng)的功率譜， 在這里，由于，a=0.95，可以算出輸出信號(hào)的方差為， 可以用留數(shù)法簡(jiǎn)單計(jì)算出來(lái)。下面對(duì)輸入輸出信號(hào)的均值、方差、相關(guān)函數(shù)及功率譜密度分別進(jìn)行討論。均值變化輸入為白噪聲，并且均值為 0，按照理論公式，可得到 下面對(duì)實(shí)際值進(jìn)行分析：輸入的隨機(jī)序列，服從標(biāo)準(zhǔn)正態(tài)分布?？梢杂孟旅娴恼Z(yǔ)句產(chǎn)生x = randn(1,500); % 產(chǎn)生題設(shè)的隨機(jī)序列,長(zhǎng)度為500 點(diǎn)系統(tǒng)的沖激響應(yīng)為，可以用下面的語(yǔ)句產(chǎn)生這個(gè)沖激信號(hào)：b=1;a=1,-0.5; % 設(shè)置濾波器的參數(shù)，b 為分子系數(shù)，a 為分母系數(shù)h=impz(b,a,20); % 得到這個(gè)系統(tǒng)的沖激響

3、應(yīng)，就是題設(shè)中的h（n）輸入信號(hào)通過(guò)線性系統(tǒng)，可以通過(guò)卷積的方法，或者用 filter 函數(shù)，y1=filter(b,a,x); % 用濾波器的方法，點(diǎn)數(shù)為500 點(diǎn)y2=conv(x,h); % 通過(guò)卷積方法得到，點(diǎn)數(shù)為519 點(diǎn)實(shí)現(xiàn)的MATLAB代碼如下：clear all;x = randn(1,500); % 產(chǎn)生題設(shè)的隨機(jī)序列,長(zhǎng)度為500點(diǎn)b=1;a=1,-0.5; % 設(shè)置濾波器的參數(shù)，b為分子系數(shù)，a為分母系數(shù)h=impz(b,a,20); % 得到這個(gè)系統(tǒng)的沖激響應(yīng)，就是題設(shè)中的h（n）y1=filter(b,a,x); % 用濾波器的方法，點(diǎn)數(shù)為500點(diǎn)subplot(2,

4、1,1); plot(y1,'r');Title('鄒先雄用濾波器的方法，點(diǎn)數(shù)為500 點(diǎn)');x = randn(1,500);y2=conv(x,h); % 通過(guò)卷積方法得到，點(diǎn)數(shù)為519點(diǎn)subplot(2,1,2); plot(y2,'b');title('鄒先雄通過(guò)卷積方法得到，點(diǎn)數(shù)為519 點(diǎn)');grid on;下面畫出兩者得到波形的區(qū)別：（為了保持一致，對(duì)y2 的輸出取前500 點(diǎn)）兩者的輸出波形近似一致，可以采用任意一個(gè)進(jìn)行分析。就采用 y1 進(jìn)行討論，輸出均值為：y1_mean=mean(y1); % 進(jìn)行時(shí)間

5、平均，求均值最終值為-0.0973，與理論的零值有一定誤差，考慮到輸入隨機(jī)序列的均值不是0，m_x=mean(x)=-0.0485,按照上面式子，得到m_y=m_xH(0)=2m_x=-0.0970理論值和實(shí)際值是非常吻合的。附運(yùn)行結(jié)果圖：*因?yàn)槭请S機(jī)序列，所以每次運(yùn)行得到y(tǒng)1和m_x的值也是隨機(jī)的，但是它們始終滿足y1=2m_x方差變化輸入信號(hào)方差的理論值就是 1，按照公式，輸出的功率譜為 下面對(duì)實(shí)際值進(jìn)行分析，用y1_var=var(y1); 求得輸出均值為 1.3598，與理論值的1.3333 有差距。如圖：自相關(guān)函數(shù)的理論與實(shí)際值理論值為： 在題設(shè)中，為白噪聲，所以 所以，輸出的自相關(guān)

6、函數(shù)理論值為 可以得到，在零點(diǎn)的值就是 1.3333，也就是輸出信號(hào)的平均功率。由MATLAb計(jì)算的結(jié)果為1.3608，這和計(jì)算結(jié)果非常接近，實(shí)際的自相關(guān)函數(shù)曲線為：clear all;x = randn(1,500); % 產(chǎn)生題設(shè)的隨機(jī)序列,長(zhǎng)度為500點(diǎn)b=1;a=1,-0.5; % 設(shè)置濾波器的參數(shù)，b為分子系數(shù)，a為分母系數(shù)h=impz(b,a,20); % 得到這個(gè)系統(tǒng)的沖激響應(yīng)，就是題設(shè)中的h（n）y1=filter(b,a,x); % 用濾波器的方法，點(diǎn)數(shù)為500點(diǎn)y2=conv(x,h); % 通過(guò)卷積方法得到，點(diǎn)數(shù)為519點(diǎn)Y3=var(y1) title('自相關(guān)

7、函數(shù)'); Ry=xcorr(y1,20,'coeff'); % 進(jìn)行歸一化的自相關(guān)函數(shù)估計(jì)，相關(guān)長(zhǎng)度為20n=-20:1:20;stem(n,Ry,'MarkerFaceColor','red');title('鄒先雄實(shí)際的自相關(guān)函數(shù)曲線');功率譜密度函數(shù)的理論與實(shí)際值對(duì)于理論的功率譜密度，可以表示為， 而對(duì)于觀測(cè)數(shù)據(jù)，可以用功率譜估計(jì)的方法得到功率譜密度。首先，采用 Welch 法估計(jì)信號(hào)的功率譜。它的原理是將數(shù)據(jù)分成等長(zhǎng)度的小段，并且允許數(shù)據(jù)的重疊，對(duì)每段進(jìn)行估計(jì)，再進(jìn)行平均，得到信號(hào)的功率譜。在Matlab 中有

8、專用函數(shù)pwelch，它的用法是：Px,f=pwelch(X,WINDOW,NOVERLAP,NFFT,Fs, 'onesided'); % window 是采用的數(shù)據(jù)窗，NOVERLAP 是重疊的數(shù)目，NFFT 是做FFT 的點(diǎn)數(shù)，F(xiàn)s 是采樣頻率，onesided 是頻率取值。針對(duì)本例，可以用下面語(yǔ)句實(shí)現(xiàn)：window=hamming(20); % 采用hanmming窗，長(zhǎng)度為20 noverlap=10; % 重疊的點(diǎn)數(shù)Nfft=512; % 做FFT的點(diǎn)數(shù)Fs=1000; % 采樣頻率，為1000Hz x = randn(1,500); % 產(chǎn)生題設(shè)的隨機(jī)序列,長(zhǎng)度為

9、500點(diǎn)b=1;a=1,-0.5; % 設(shè)置濾波器的參數(shù)，b為分子系數(shù)，a為分母系數(shù)h=impz(b,a,20); % 得到這個(gè)系統(tǒng)的沖激響應(yīng)，就是題設(shè)中的h（n）y1=filter(b,a,x);y1_mean=mean(y1); % 進(jìn)行時(shí)間平均，求均值y1_var=var(y1); % 進(jìn)行時(shí)間平均，求方差Ry=xcorr(y1,20,'coeff'); % 進(jìn)行歸一化的自相關(guān)函數(shù)估計(jì)，相關(guān)長(zhǎng)度為20Py,f=pwelch(y1,window,noverlap,Nfft,Fs, 'onesided'); % 估計(jì)功率譜密度f(wàn)=-fliplr(f) f(1:

10、end); % 構(gòu)造一個(gè)對(duì)稱的頻率，范圍是-Fs/2, Fs/2 Py=-fliplr(Py) Py(1:end); % 對(duì)稱的功率譜plot(f,10*log10(abs(Py),'r');title('鄒先雄實(shí)際功率譜密度曲線');grid on;最后，得到的估計(jì)值為根據(jù)上述值，可以計(jì)算出理論的功率,由于 可以用下面的語(yǔ)句實(shí)現(xiàn)：w=2*pi*f/Fs; % 轉(zhuǎn)化到數(shù)字域上面H=(1+0.25-2*0.5*cos(w).(-1);% 系統(tǒng)函數(shù)模平方Gy=H/max(H); % 歸一化處理Gy=10*log10(Gy); % 化成dB形式plot(f,Gy,&#

11、39;b');title('鄒先雄理論功率譜密度曲線');grid on;畫出的圖形見(jiàn)下圖：這是理論的功率譜密度。為了方便顯示，將兩幅圖畫在一起，便于比較。window=hamming(20); % 采用hanmming窗，長(zhǎng)度為20 noverlap=10; % 重疊的點(diǎn)數(shù)Nfft=512; % 做FFT的點(diǎn)數(shù)Fs=1000; % 采樣頻率，為1000Hz x = randn(1,500); % 產(chǎn)生題設(shè)的隨機(jī)序列,長(zhǎng)度為500點(diǎn)b=1;a=1,-0.5; % 設(shè)置濾波器的參數(shù)，b為分子系數(shù)，a為分母系數(shù)h=impz(b,a,20); % 得到這個(gè)系統(tǒng)的沖激響應(yīng)，就是

12、題設(shè)中的h（n）y1=filter(b,a,x);y1_mean=mean(y1); % 進(jìn)行時(shí)間平均，求均值y1_var=var(y1); % 進(jìn)行時(shí)間平均，求方差Ry=xcorr(y1,20,'coeff'); % 進(jìn)行歸一化的自相關(guān)函數(shù)估計(jì)，相關(guān)長(zhǎng)度為20Py,f=pwelch(y1,window,noverlap,Nfft,Fs, 'onesided'); % 估計(jì)功率譜密度f(wàn)=-fliplr(f) f(1:end); % 構(gòu)造一個(gè)對(duì)稱的頻率，范圍是-Fs/2, Fs/2 Py=-fliplr(Py) Py(1:end); % 對(duì)稱的功率譜Py=Py/m

13、ax(Py); % 歸一化處理w=2*pi*f/Fs; % 轉(zhuǎn)化到數(shù)字域上面H=(1+0.25-2*0.5*cos(w).(-1);% 系統(tǒng)函數(shù)模平方Gy=H/max(H); % 歸一化處理Gy=10*log10(Gy); % 化成dB形式plot(f,10*log10(abs(Py),'r',f,Gy,'b');title('鄒先雄實(shí)際功率譜和理論功率譜擬合');legend(' ','實(shí)際值','理論值');grid on;結(jié)果為：從結(jié)果上可以看出來(lái)，兩者存在著比較大的差距，這是由于輸入隨機(jī)序列

14、的功率譜并不是常數(shù)的緣故，也就是輸入不是嚴(yán)格的白噪聲，所以會(huì)出現(xiàn)波動(dòng)。當(dāng)隨著數(shù)據(jù)值的增加，擬合的程度會(huì)有所改善。3、實(shí)驗(yàn)內(nèi)容假定一具有單位方差的抽樣序列X(n)的白噪聲隨機(jī)過(guò)程X(t)通過(guò)一脈沖響應(yīng)為 的線性濾波器，利用 matlab 工具繪出輸入輸出信號(hào)的均值、方差、相關(guān)函數(shù)及功率譜密度。實(shí)現(xiàn)的MATLAB代碼和結(jié)果如下:clear all;x = randn(1,500); % 產(chǎn)生題設(shè)的隨機(jī)序列,長(zhǎng)度為500點(diǎn)b=1;a=1,-0.6; % 設(shè)置濾波器的參數(shù)，b為分子系數(shù)，a為分母系數(shù)h=impz(b,a,20); % 得到這個(gè)系統(tǒng)的沖激響應(yīng)，就是題設(shè)中的h（n）y1=filter(b,

15、a,x); % 用濾波器的方法，點(diǎn)數(shù)為500點(diǎn)figure(1)subplot(2,1,1); plot(y1,'r');title('鄒先雄用濾波器的方法，點(diǎn)數(shù)為500 點(diǎn)');y2=conv(x,h); % 通過(guò)卷積方法得到，點(diǎn)數(shù)為519點(diǎn)subplot(2,1,2); plot(y2,'b');title('鄒先雄通過(guò)卷積方法得到，點(diǎn)數(shù)為519 點(diǎn)');Y2=mean(y1) %進(jìn)行時(shí)間平均，求均值，為理論值 m_x=mean(x) %輸出的實(shí)際值可以通過(guò)2m_x計(jì)算Y3=var(y1) figure(2) Ry=xcorr

16、(y1,20,'coeff'); % 進(jìn)行歸一化的自相關(guān)函數(shù)估計(jì)，相關(guān)長(zhǎng)度為20n=-20:1:20;stem(n,Ry,'MarkerFaceColor','red');title('鄒先雄實(shí)際的自相關(guān)函數(shù)曲線');%實(shí)際功率譜密度window=hamming(20); % 采用hanmming窗，長(zhǎng)度為20 noverlap=10; % 重疊的點(diǎn)數(shù)Nfft=512; % 做FFT的點(diǎn)數(shù)Fs=1000; % 采樣頻率，為1000Hz x = randn(1,500); % 產(chǎn)生題設(shè)的隨機(jī)序列,長(zhǎng)度為500點(diǎn)b=1;a=1,-0.

17、5; % 設(shè)置濾波器的參數(shù)，b為分子系數(shù)，a為分母系數(shù)h=impz(b,a,20); % 得到這個(gè)系統(tǒng)的沖激響應(yīng)，就是題設(shè)中的h（n）y1=filter(b,a,x);y1_mean=mean(y1); % 進(jìn)行時(shí)間平均，求均值y1_var=var(y1); % 進(jìn)行時(shí)間平均，求均值Ry=xcorr(y1,20,'coeff'); % 進(jìn)行歸一化的自相關(guān)函數(shù)估計(jì)，相關(guān)長(zhǎng)度為20Py,f=pwelch(y1,window,noverlap,Nfft,Fs, 'onesided'); % 估計(jì)功率譜密度f(wàn)=-fliplr(f) f(1:end); % 構(gòu)造一個(gè)對(duì)稱的

18、頻率，范圍是-Fs/2, Fs/2 Py=-fliplr(Py) Py(1:end); % 對(duì)稱的功率譜figure(3)plot(f,10*log10(abs(Py),'r');title('鄒先雄實(shí)際功率譜密度');grid on;%理論功率譜密度w=2*pi*f/Fs; % 轉(zhuǎn)化到數(shù)字域上面H=(1+0.25-2*0.5*cos(w).(-1);% 系統(tǒng)函數(shù)模平方Gy=H/max(H); % 歸一化處理Gy=10*log10(Gy); % 化成dB形式figure(4)plot(f,Gy,'b');title('鄒先雄理論功率譜密度');grid on;下圖為隨機(jī)產(chǎn)生的波形，以此來(lái)求相應(yīng)的均值、方差、相關(guān)函數(shù)及功率譜密度。自相關(guān)函數(shù)波形功率譜密度波形如下圖截的部分程序及運(yùn)行結(jié)果，其中y2為輸出的均值，y3為方差4.實(shí)驗(yàn)總結(jié)： 隨機(jī)信號(hào)實(shí)驗(yàn)綜合性強(qiáng)，運(yùn)用到了隨機(jī)信號(hào)分析，數(shù)字信號(hào)處理，概率論，matlab等課程知識(shí)，讓我們嘗試綜合運(yùn)用理論知識(shí)，這樣，我們能夠更深刻的理解理論知