如何在Matlab/Simulink中生成指定的白噪声和有色噪声

在进行系统仿真时，经常需要利用白噪声和有色噪声作为系统输入，在Matlab和Simulink环境下提供了多种方式生成白噪声或者有色噪声。在这里对自己常用的一些方法进行简单总结。欢迎一起交流并提供更多的思路。

白噪声和有色噪声基础知识
在Matlab/Simulink中生成单位功率谱密度的白噪声
在Matlab/Simulink中利用成形滤波器生成有色噪声
功率谱密度分析

1. 白噪声和有色噪声基础知识

在学术声，白噪声和有色噪声的定义如下：

白噪声（white noise）是指功率谱密度在整个频域内是常数的噪声。所有频率具有相同能量密度的随机噪声称为白噪声。

有色噪声（ coloured noise）是指功率谱密度函数不平坦的噪声。大多数的噪声的频谱主要都是非白色频谱，通过信道的白噪声受信道频率的影响而变为有色的。

在物理世界中，任何系统都是有限带宽的。而根据白噪声的定义，其带宽是无穷大的，意味着其能量是无限的，这显然是不现实的。所以，在数学分析或者仿真时，我们常认为在有限带宽内具有平坦功率谱密度的信号就是白噪声。

$E[y(t)]=E[x(t)]\int_0^{\infty}h(\lambda)d\lambda$
式中， $h(\lambda)(\lambda>0)$ 为 $G(s)$ 的脉冲响应函数，则输入信号 $x(t)$ 和输出信号 $y(t)$ 的功率谱密度具有如下对应关系：
$S_y(\omega)=|G(j\omega)|^2S_x(\omega)$
式中， $G(j\omega)$ 为传递函数 $G(s)$ 的频率响应。因此，可以通过将白噪声信号 $x(t)$ 通过线性传递函数 $G(s)$ ，则可以得到有色噪声信号 $y(t)$ ，如下所示：
$S_y(\omega)=|G(j\omega)|^2S_x(\omega)=|G(j\omega)|^2 \cdot 1=G^*(j\omega)G(j\omega)$
$G^*(j\omega)$ 为 $G(j\omega)$ 的复共轭。 $G(s)$ 就是成形滤波器的传递函数形式。

下图给出了白噪声和有色噪声的功率谱密度对比图：

下面结合实际代码，介绍如何在Matlab/Simulink中分别生成指定功率谱密度的白噪声和有色噪声。

2. 在Matlab/Simulink中生成单位功率谱密度的白噪声

功率谱密度与频谱图一样，关于 $f=0$ 对称，即 $S(f)=S(-f)$ 。因此对于频率 $f\in (-\infty, +\infty)$ 的功率谱密度，称之为双边功率谱密度（double-side power spectrum density）；对于频率 $f\in(0,+\infty)$ 则称为单边功率谱密度密度（single-side power spectrum density）。由双边功率谱密度转换为单边功率谱密度时需要乘以系数2:
$\int_{-\infty}^{+\infty}S_{xx}(f)df=2\int_0^{+\infty}S_{x}(f)df$
下面分别介绍两种功率谱密度在Matlab/Simulink中的生成方式。

2.1 双边功率谱密度为1的白噪声

在Matlab中可以用随机数函数生成白噪声，例如randn()函数。根据白噪声的定义，均值为0，需要确定的只剩随机序列的方差。在前述理论基础中已经介绍了，在工程上，我们假设在有限带宽上具有常值功率谱密度的信号就是白噪声，因此有：
$\begin{aligned} \sigma^2_x&=2\int^F_0S_{xx}df=2F\cdot \frac{1}{F}\int^F_0S_{xx}(f)df=2F\cdot E[S_{xx}(f)]\\ &=\int^F_0S_{x}df=F\cdot \frac{1}{F}\int^F_0S_{x}(f)df=F\cdot E[S_{x}(f)] \end{aligned}$
式中， $S_{xx}(f)$ 表示双边功率谱密度； $S_x(f)$ 表示单边功率谱密度； $F$ 表示系统带宽。

   若系统的采样频率$f_s=10\text{ Hz}$，根据奈奎斯特采样定理，系统最高频率或带宽则为$F = f_s/2=5 \text{ Hz}$。所以对于双边功率谱密度为$S_{xx}(f)=1$的信号，对应的噪声序列的方差为$\sigma^2 = 2\times5\times1=10$。对应的Matlab代码如下：

noise_power=1;
fs=10;
ts=1/fs;
t1_m = 0:1/fs:100;
x1_m = sqrt(2*fs/2*noise_power)*randn(length(t1_m),1); %fs/2 is Nyquist frequency

在Simulink中则提供了更加方便的方法。采用有限带宽白噪声模块（Band-limited White Noise）可以很方便的生成指定功率谱密度的白噪声。该模块需要设置的参数和定义如下：

Noise power：白噪声功率谱密度的幅值
Sample time：采样时间
Seed：随机数种子

需要注意的是该模块的中的功率谱密度的定义是双边功率谱密度。因此可以在该模块中赋如下值：

两种方法生成的生成的白噪声的仿真结果对比如下，可以看出所得到的结果是一致的。

利用mean()函数计算功率谱密度的平均值，得到的结果分别为： $S_{xx-simulink}=1.0161$ 和 $S_{xx-matlab}=0.9850$ 。可见仿真结果与预期是一致的。

2.2 单边功率谱密度为1的白噪声

有了前述分析，单边功率谱密度为1的白噪声就可以类似得到了。

首先，在Matlab中，唯一要确定的就是噪声序列的方差。同样，若系统的采样频率 $f_s=10\text{ Hz}$ ，根据奈奎斯特采样定理，系统最高频率或带宽则为 $F = f_s/2=5 \text{ Hz}$ 。所以对于单边功率谱密度为 $S_{x}(f)=1$ 的信号，对应的噪声序列的方差为 $\sigma^2 = F\cdot S_x(f)=5$ 。对应的Matlab代码如下：

noise_power=1;
fs=10;
ts=1/fs;
t1_m = 0:1/fs:100;
x1_m = sqrt(fs/2*(noise_power))*randn(length(t1_m),1); %fs/2 is Nyquist frequency

在Simulink中，需要注意有限带宽白噪声模块的noise power值。由于在Simulink中，该值的定义为双边功率谱密度，根据双边功率谱密度和单边功率谱密度的关系，noise power需要设定为0.5，如下所示：

两种方法生成的生成的白噪声的仿真结果对比如下，可以看出所得到的结果是一致的。

利用mean()函数计算功率谱密度的平均值，得到的结果分别为： $S_{x-simulink}=1.0141$ 和 $S_{x-matlab}=0.9902$ 。可见仿真结果与预期是一致的。

3. 在Matlab/Simulink中利用成形滤波器生成有色噪声

根据基础知识，将单位功率谱密度的白噪声通过成形滤波器可以得到有色噪声。成型滤波器可以通过传递函数建模得到，白噪声的建立已经分析过了，因此剩余的工作就十分简单了。

利用2.1中生成的双边功率谱密度为1的白噪声，在Matlab中定义传递函数如下所示的成形滤波器。
$G(s)=\frac{(s+1)^2}{(100s+1)^2}$
可以看出该成形滤波器也可以看作是一个低通滤波器。再利用lsim()函数则可以仿真得到有色噪声了，代码如下：

noise_power=1;
fs=10;
ts=1/fs;

t1_m = 0:1/fs:10000;
x1_m = sqrt(2*fs/2*noise_power)*randn(length(t1_m),1); %fs/2 is Nyquist frequency

s = tf('s');
shapeFilter = (1e0*s+1)^2/(1e2*s+1)^2;
y = lsim(shapeFilter,x1_m,t1_m);

在simulink中可以更加直观的进行仿真。搭建如下模块，则可以进行仿真得到有色噪声：

有色噪声仿真结构-水印.png

对比仿真结果如下所示。可以看出，两种方法得到的结果一致。仿真得到的有色噪声也与成形滤波器符合，证明了仿真的准确性。同时，中间的图也表明输入白噪声的高频变化被成形滤波器的低通滤波特性所滤去了。

4. 功率谱密度分析

功率谱密度是随机过程中无法绕开的概念。关于功率谱密度的相关介绍可以参考这篇文章--如何理解随机振动的功率谱密度？。这篇文章非常形象生动，我就是从这篇文章开始理解了什么是功率谱密度，在此也表示对作者深深地感谢。

在 Matlab中（simulink中也有相应模块，可以自行搜索power spectrum density）提供了三种计算功率谱密度的方法：

利用快速傅立叶变换，在帮助文件中搜索Power Spectral Density Estimates Using FFT，可以学习相关算法
pwelch()函数
periodogram()函数

本次所有的功率谱密度计算均是采用periodogram()函数。下一篇文章计划结合实际例子和代码，总结功率谱密度的相关概念和知识，欢迎关注。