随机信号分析实验
一、实验目的：
⑴ 了解随机信号自身的特性，包括均值（数学期望）、均方值、方差、相关函数、概率密度、频谱及功率谱密度等。
（2）了解色噪声的基本概念和分析方法，掌握用matlab、c\\c++软件仿真和分析色噪声的方法。
⑶ 掌握随机信号的分析方法。
二、实验原理：
我们把除了白噪声之外的所有噪声都称为有色噪声。就像白光一样，除了白光就是有色光。
色噪声中有几个典型：
⑴ 粉红噪声。粉红噪音是自然界最常见的噪音，简单说来，粉红噪音的频率分量功率主要分布在中低频段。从波形角度看，粉红噪音是分形的，在一定的范围内音频数据具有相同或类似的能量。从功率（能量）的角度来看，粉红噪音的能量从低频向高频不断衰减，曲线为1/f，通常为每8度下降3分贝。粉红噪声的能量分布在任一同比例带宽中是相等的！比如常见的三分之一倍程频带宽100Hz的范围.2__112和1000Hz的2__1120是相等的。
在给定频率范围内(不包含直流成分)，随着频率的增加，其功率密度每倍频程下降3dB(密度与频率成反比)。每倍频的功率相同，但要产生每倍频程3dB的衰减非常困难,因此，没有纹波的粉红噪声在现实中很难找到。
粉红噪声低频能下降到接近0Hz(不包括0Hz)高频端能上到二十几千赫，而且它在等比例带宽内的能量是相等的(误差只不过0.1dB左右)。 粉红噪声的功率普密度图：
⑵ 红噪声(海洋学概念)。这是有关海洋环境的一种噪声，由于它是有选择地吸收较高的频率，因此称之为红噪声。
⑶ 橙色噪声。该类噪声是准静态噪声，在整个连续频谱范围内，功率谱有限且零功率窄带信号数量也有限。这些零功率的窄带信号集中于任意相关音符系统的音符频率中心上。由于消除了所有的合音，这些剩余频谱就称为橙色音符。 ⑷ 蓝噪声。在有限频率范围内，功率密度随频率的增加每倍频增长3dB(密
度正比于频率)。对于高频信号来说，它属于良性噪声。 ⑸ 紫噪声。在有限频率范围内，功率密度随频率的增加每倍频增长6dB(密度正比于频率的平方值)。 ⑹ 灰色噪声。该噪声在给定频率范围内，类似于心理声学上的等响度曲线(如反向的A-加权曲线)，因此在所有频率点的噪声电平相同。 ⑺ 棕色噪声。在不包含直流成分的有限频率范围内，功率密度随频率的增加每倍频下降6dB(密度与频率的平方成反比)。该噪声实际上是布朗运动产生的噪声，它也称为随机飘移噪声或醉鬼噪声。 ⑻ 黑噪声(静止噪声)包括：
① 有源噪声控制系统在消除了一个现有噪声后的输出信号。 ② 在20kHz以上的有限频率范围内，功率密度为常数的噪声，一定程度上它类似于超声波白噪声。 这种黑噪声就像“黑光”一样，由于频率太高而使人们无法感知，但它对你和你周围的环境仍然有影响。
三、实验内容：
⑴ 用matlab或c/c++语言编写和仿真程序。
⑵ 产生粉红色噪声和高斯色噪声：让高斯白噪声通过低通、带通、高通滤波器中的任意一个就可以产生高斯色噪声。让高斯白噪声通过每倍频程衰减3dB的衰减滤波器的滤波器就可以产生粉红噪声。
⑶ 对粉红色噪声和高斯色噪声进行相关分析和谱分析。计算粉红色噪声、高斯色噪声的均值、均方值、方差、相关函数、概率密度、频谱及功率谱密度、相关函数。
⑷ 所有结果均用图示法来表示，能读出具体值。
四、实验中产生的各种波形：
白噪声及其各种波
白噪声是随机的，所以它的均值为零，频谱在所有频率上都有分量。功率
谱密度近似均匀分布。因为它不具有相关性，所以自相关函数趋近于零。白噪声有波动所以方差不为零。
低通滤波器
试验中采用白噪声通过低通滤波器的方法来产生高斯色噪声，滤波器通带截止频率为30Hz，阻带截止频率为40Hz。
高斯色噪声及其各种波
色噪声的功率谱不是均匀分布的，它的均值不为零，具有一定的波动所以
方差也不为零。它的前后时刻是相关的，所以自相关函数不为零。
倍频程衰减滤波器
倍频程衰减3dB滤波器的设计参考《DSP generation of Pink (1/f) Noise》 作者Robert Bristow-Johnson。
粉红噪声及其各种波
色噪声的功率谱不是均匀分布的，它的均值不为零，具有一定的波动所以方差也不为零。它的前后时刻是相关的，所以自相关函数不为零。
五、附：实验参考书目及Matlab程序：
参考书目：
（1）楼顺天编著.基于MATLAB的系统分析与设计. 西安电子科技大学出版社,2002 （2）[美]Sanjit KMitra编著.Digital Signal Processing——A Computer-Based Approach
电子工业出版社，2010
（3）苏晓生编著.掌握Matlab6.0及其工程应用. 科学出版社，2004
（4）徐小兵 沈勇 邬宁.IIR 数字粉红噪声滤波器的优化设计。电声技术，2005-12
（5）Robert Bristow-Johnson.《DSP generation of Pink (1/f) Noise》
试验程序如下：
%产生高斯白噪声 clf;
y=wgn(1,1024,0); t=0:1023;
y1=mean(y);%均值01 y2=var(y);%方差
y3=y2+y1.*y1;%均方值
[y4,lag]=xcorr(y,'unbiased');%自相关函数 [f1,y5] = ksdensity(y);%概率密度 f=(0:length(y)-1)'/length(y)*1024; Y=fft(y);
y6=abs(Y);%频谱
y7=Y.*conj(Y)/1024;%功率谱密度 figure(1);
subplot(2,4,1);plot(t,y);
title('高斯白噪声');axis([0 1024 -5 5]); subplot(2,4,2);plot(t,y1);
title('高斯白噪声均值');axis([0 1024 -2 2]); subplot(2,4,3);plot(t,y2);
title('高斯白噪声方差');axis([0 1024 -2 2]); subplot(2,4,4);plot(t,y3);
title('高斯白噪声均方值');axis([0 1024 -2 2]); subplot(2,4,5);plot(lag,y4);
title('高斯白噪声自相关函数');axis([-1024 1024 -1 1]); subplot(2,4,6);plot(y5,f1); title('概率密度');
subplot(2,4,7);plot(f,y6);
title('高斯白噪声频谱');axis([0 1024 0 80]);
subplot(2,4,8);plot(f,y7);
title('高斯白噪声功率谱密度');axis([0 1024 0 8]); %低通滤波器
Wp=2*pi*30;Ws=2*pi*40;Rp=0.5;Rs=40;fs=100;W=2*pi*fs; [N,Wn]=buttord(2*Wp/W,2*Ws/W,Rp,Rs); [b,a]=butter(N,Wn);
[h,f]=freqz(b,a,1000,fs);
figure(2);plot(f,abs(h)); xlabel('f/Hz');ylabel('|H(jf)|'); axis([0 100 0 1.2]);grid on;
title('低通滤波器幅频响应');
%生成高斯色噪声
gss=filter(b,a,y);%滤波产生高斯色噪声 gss1=mean(gss);%均值 gss2=var(gss);%方差
gss3=gss2+gss1.*gss1;%均方值
[gss4,lag]=xcorr(gss,'unbiased');%自相关函数 [f1,gss5]=ksdensity(gss);%概率密度
f=(0:length(gss)-1)/length(gss)*1024; GSS=fft(gss);
gss6=abs(GSS);%频谱
gss7=GSS.*conj(GSS)/1024;%功率谱密度 figure(3);
subplot(2,4,1);plot(t,gss);
title('高斯色噪声');axis([0 1024 -5 5]); subplot(2,4,2);plot(t,gss1);
title('高斯色噪声均值');axis([0 1024 -1 1]); subplot(2,4,3);plot(t,gss2);
title('高斯色噪声方差');axis([0 1024 -0.5 1.5]); subplot(2,4,4);plot(t,gss3);
title('高斯色噪声均方值');axis([0 1024 -0.5 1.5]); subplot(2,4,5);plot(lag,gss4);
title('高斯色噪声自相关函数');axis([-1024 1024 -0.5 1]); subplot(2,4,6);plot(gss5,f1); title('高斯色噪声概率密度'); subplot(2,4,7);plot(f,gss6);
title('高斯色噪声频谱');axis([0 1024 0 100]); subplot(2,4,8);plot(f,gss7);
title('高斯色噪声功率谱密度');axis([0 1024 0 8]);
%倍频程衰减3dB滤波器
%倍频程衰减3dB滤波器的设计参考
%《DSP generation of Pink (1/f) Noise》 %作者Robert Bristow-Johnson
a=[0.98443604 0.83392334 0.07568359]; b=[0.99572754 0.947909 0.53567505]; Hz = zpk(a,b,1,1/44100); Hz2=get(tf(Hz))
num = cell2mat(Hz2.num()); den = cell2mat(Hz2.den(1)); [H F]=freqs(den,num);
figure(4);semilogx(F,20*log(abs(H)));
xlabel('F/Hz');ylabel('|H(jf)|/dB');grid on; title('倍频程衰减滤波器幅频响应');
%产生粉红噪声
pn=filter(num,den,y);%滤波产生粉红色噪声 pn1=mean(pn);%粉红噪声均值 pn2=var(pn);%粉红噪声方差
pn3=pn2+pn1.*pn1;%粉红噪声均方值
[pn4,lag]=xcorr(pn,'unbiased');%粉红噪声自相关函数 [f1,pn5]=ksdensity(pn);%粉红噪声概率密度 f=(0:length(pn)-1)/length(pn)*1024; PN=fft(pn);
pn6=abs(PN);%粉红噪声频谱
pn7=PN.*conj(PN)/1024;%粉红噪声功率谱密度 figure(5);
subplot(2,4,1);plot(t,pn);
title('粉红噪声');axis([0 1024 -8 8]); subplot(2,4,2);plot(t,pn1);
title('粉红噪声均值');axis([0 1024 -1 1]); subplot(2,4,3);plot(t,pn2);
title('粉红噪声方差');axis([0 1024 0 5]); subplot(2,4,4);plot(t,pn3);
title('粉红噪声均方值');axis([0 1024 0 5]); subplot(2,4,5);plot(lag,pn4);
title('粉红噪声自相关函数');axis([-1024 1024 -1 2]); subplot(2,4,6);plot(pn5,f1); title('粉红噪声概率密度');
subplot(2,4,7);plot(f,pn6);
title('粉红噪声频谱');axis([-50 1080 0 350]); subplot(2,4,8);plot(f,pn7);
title('粉红噪声功率谱密度');axis([-50 1080 0 110]);