白噪聲通過LTI-的仿真

1、實驗2 白噪聲通過LTI 的仿真1、實驗目的了解白噪聲通過LTI 系統(tǒng)的原理與處理方法，學會運用Matlab 函數對隨機過程進行均值、相關函數和功率譜的估計,并且通過實驗分析理論分析與實驗結果之間的差異。2、實驗原理假定一具有單位方差的抽樣序列X(n)的白噪聲隨機過程X(t)通過一脈沖響應為 的線性濾波器，繪出輸入輸出信號的均值、方差、相關函數及功率譜密度。設系統(tǒng)沖激響應為h(n)，傳遞函數 ，或者用Z 變換，結果為 。輸入為 X (n)，輸出為 ，均值關系： ，若平穩(wěn)有， 自相關函數關系， ，當是平穩(wěn)時候，有 題目中假設為白噪聲，可以根據白噪聲的性質進行理論計算。白噪聲的自相關函數， 這里，

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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