色噪声的产生及MATLAB实现

色噪声

色噪声是指功率谱不平坦的噪声，与之相对应的是白噪声，白噪声定义为在无限频率范围内功率谱密度为常数的信号。一种常见的色噪声产生方法就是利用白噪声通过滤波器，使其功率谱在对应的范围内不再保持常数。

ARMA（自回归滑动平均）模型

如果系数 a 1 , a 2 , ⋯ , a p {{a}_{1}},{{a}_{2}},\cdots ,{{a}_{p}} a1,a2,⋯,ap和 b 1 , b 2 , ⋯ , b q {{b}_{1}},{{b}_{2}},\cdots ,{{b}_{q}} b1,b2,⋯,bq不全为零，并且 a p ≠ 0 {{a}_{p}}\ne 0 ap=0， b q ≠ 0 {{b}_{q}}\ne 0 bq=0，那么此时系统既存在零点也存在极点。当序列 w ( n ) w\left( n \right) w(n)通过上述系统后，那么系统的输出 y ( n ) y\left( n \right) y(n)与输入 w ( n ) w\left( n \right) w(n)之间的关系可以用下列查分方程表示
y ( n ) = − a 1 y ( n − 1 ) − a 2 y ( n − 2 ) − ⋯ − a p y ( n − p ) + w ( n ) + b 1 w ( n − 1 ) + b 2 w ( n − 2 ) + ⋯ + b q w ( n − q ) = ∑ k = 0 q b k w ( n − k ) − ∑ k = 1 p a k y ( n − k ) \begin{aligned} & y\left( n \right)=-{{a}_{1}}y\left( n-1 \right)-{{a}_{2}}y\left( n-2 \right)-\cdots -{{a}_{p}}y\left( n-p \right) \\ & \ \ +w\left( n \right)+{{b}_{1}}w\left( n-1 \right)+{{b}_{2}}w\left( n-2 \right)+\cdots +{{b}_{q}}w\left( n-q \right) \\ & =\sum\limits_{k=0}^{q}{{{b}_{k}}w\left( n-k \right)}-\sum\limits_{k=1}^{p}{{{a}_{k}}y\left( n-k \right)} \end{aligned} y(n)=−a1y(n−1)−a2y(n−2)−⋯−apy(n−p) +w(n)+b1w(n−1)+b2w(n−2)+⋯+bqw(n−q)=k=0∑qbkw(n−k)−k=1∑paky(n−k)

系统输出功率谱密度函数

S y ( Ω ) = S w ( Ω ) ∣ ∑ k = 0 q b k e − j Ω k ∣ 2 ∣ ∑ k = 0 p a k e − j Ω k ∣ 2 {{S}_{y}}\left( \Omega \right)={{S}_{w}}\left( \Omega \right)\frac{{{\left| \sum\limits_{k=0}^{q}{{{b}_{k}}{{e}^{-j\Omega k}}} \right|}^{2}}}{{{\left| \sum\limits_{k=0}^{p}{{{a}_{k}}{{e}^{-j\Omega k}}} \right|}^{2}}} Sy(Ω)=Sw(Ω)∣∣∣∣k=0∑pake−jΩk∣∣∣∣2∣∣∣∣k=0∑qbke−jΩk∣∣∣∣2

系统输出自相关函数

R y ( m ) = R y ( n 1 − n 2 ) = E { y ( n 2 ) [ ∑ k = 0 q b k w ( n 1 − k ) − ∑ k = 1 p a k y ( n 1 − k ) ] } = ∑ k = 0 q b k R w y ( m − k ) − ∑ k = 1 p a k R y ( m − k ) \begin{aligned} & {{R}_{y}}\left( m \right)={{R}_{y}}\left( {{n}_{1}}-{{n}_{2}} \right)=E\left\{ y\left( {{n}_{2}} \right)\left[ \sum\limits_{k=0}^{q}{{{b}_{k}}w\left( {{n}_{1}}-k \right)-\sum\limits_{k=1}^{p}{{{a}_{k}}y\left( {{n}_{1}}-k \right)}} \right] \right\} \\ & =\sum\limits_{k=0}^{q}{{{b}_{k}}{{R}_{wy}}\left( m-k \right)-\sum\limits_{k=1}^{p}{{{a}_{k}}{{R}_{y}}\left( m-k \right)}} \end{aligned} Ry(m)=Ry(n1−n2)=E{y(n2)[k=0∑qbkw(n1−k)−k=1∑paky(n1−k)]}=k=0∑qbkRwy(m−k)−k=1∑pakRy(m−k)
如果输入是白噪声序列的话，那么系统的输出的功率谱和自相关函数可以表示为
S y ( Ω ) = σ w 2 ∣ ∑ k = 0 q b k e − j Ω k ∣ 2 ∣ ∑ k = 0 p a k e − j Ω k ∣ 2 {{S}_{y}}\left( \Omega \right)=\sigma _{w}^{2}\frac{{{\left| \sum\limits_{k=0}^{q}{{{b}_{k}}{{e}^{-j\Omega k}}} \right|}^{2}}}{{{\left| \sum\limits_{k=0}^{p}{{{a}_{k}}{{e}^{-j\Omega k}}} \right|}^{2}}} Sy(Ω)=σw2∣∣∣∣k=0∑pake−jΩk∣∣∣∣2∣∣∣∣k=0∑qbke−jΩk∣∣∣∣2
R y ( m ) = R y ( n 1 − n 2 ) = E { y ( n 2 ) [ ∑ k = 0 q b k w ( n 1 − k ) − ∑ k = 1 p a k y ( n 1 − k ) ] } = ∑ k = 0 q b k R w y ( m − k ) − ∑ k = 1 p a k R y ( m − k ) = ∑ k = 0 q b k σ w 2 h ( m − k ) − ∑ k = 1 p a k R y ( m − k ) \begin{aligned} & {{R}_{y}}\left( m \right)={{R}_{y}}\left( {{n}_{1}}-{{n}_{2}} \right)=E\left\{ y\left( {{n}_{2}} \right)\left[ \sum\limits_{k=0}^{q}{{{b}_{k}}w\left( {{n}_{1}}-k \right)-\sum\limits_{k=1}^{p}{{{a}_{k}}y\left( {{n}_{1}}-k \right)}} \right] \right\} \\ & =\sum\limits_{k=0}^{q}{{{b}_{k}}{{R}_{wy}}\left( m-k \right)-\sum\limits_{k=1}^{p}{{{a}_{k}}{{R}_{y}}\left( m-k \right)}} \\ & \text{=}\sum\limits_{k=0}^{q}{{{b}_{k}}\sigma _{w}^{2}h\left( m-k \right)-\sum\limits_{k=1}^{p}{{{a}_{k}}{{R}_{y}}\left( m-k \right)}} \end{aligned} Ry(m)=Ry(n1−n2)=E{y(n2)[k=0∑qbkw(n1−k)−k=1∑paky(n1−k)]}=k=0∑qbkRwy(m−k)−k=1∑pakRy(m−k)=k=0∑qbkσw2h(m−k)−k=1∑pakRy(m−k)
其中 h ( n ) h\left( n \right) h(n)为系统的冲激响应，其表达式为
h ( n ) = { − ∑ k = 1 p a k h ( n − k ) + b n , 0 ≤ n ≤ q − ∑ k = 1 p a k h ( n − k ) , n > q 0 , n < 0 h\left( n \right)=\left\{ \begin{aligned} & -\sum\limits_{k=1}^{p}{{{a}_{k}}h\left( n-k \right)}+{{b}_{n}},\ \ \ \ \ \ \ 0\le n\le q \\ & -\sum\limits_{k=1}^{p}{{{a}_{k}}h\left( n-k \right),\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ n>q} \\ & \ \ \ \ \ \ \ \ \ \ \ 0,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ n<0 \\ \end{aligned} \right. h(n)=⎩⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎧−k=1∑pakh(n−k)+bn, 0≤n≤q−k=1∑pakh(n−k), n>q 0, n<0

AR模型

当系数 b 1 , b 2 , ⋯ , b q {{b}_{1}},{{b}_{2}},\cdots ,{{b}_{q}} b1,b2,⋯,bq 为0且 a p ≠ 0 {{a}_{p}}\ne 0 ap=0 时，系统为全极点系统，
系统的输出可以表示为
y ( n ) = w ( n ) − ∑ k = 1 p a k y ( n − k ) y\left( n \right)=w\left( n \right)-\sum\limits_{k=1}^{p}{{{a}_{k}}y\left( n-k \right)} y(n)=w(n)−k=1∑paky(n−k)
那么系统的传输函数和冲激响应可以表示为
H ( j Ω ) = 1 1 + ∑ k = 1 p a k e − j Ω k = 1 ∑ k = 0 p a k e − j Ω k H\left( j\Omega \right)=\frac{1}{1+\sum\limits_{k=1}^{p}{{{a}_{k}}{{e}^{-j\Omega k}}}}=\frac{1}{\sum\limits_{k=0}^{p}{{{a}_{k}}{{e}^{-j\Omega k}}}} H(jΩ)=1+k=1∑pake−jΩk1=k=0∑pake−jΩk1
h ( n ) = { − ∑ k = 1 p a k h ( n − k ) + δ ( n ) , n ≥ 0 0 , n < 0 h\left( n \right)=\left\{ \begin{aligned} & -\sum\limits_{k=1}^{p}{{{a}_{k}}h\left( n-k \right)+\delta \left( n \right),\ \ \ \ n\ge 0} \\ & 0,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ n<0 \\ \end{aligned} \right. h(n)=⎩⎪⎪⎨⎪⎪⎧−k=1∑pakh(n−k)+δ(n), n≥00, n<0
由上面的表达式可以看出，AR（p）模型是AR（p,q）模型经过q=0退化来的，那么系统输出的功率谱密度函数和自相关函数可以表示为
S y ( Ω ) = σ w 2 ∣ ∑ k = 0 p a k e − j Ω k ∣ 2 {{S}_{y}}\left( \Omega \right)=\frac{\sigma _{w}^{2}}{{{\left| \sum\limits_{k=0}^{p}{{{a}_{k}}{{e}^{-j\Omega k}}} \right|}^{2}}} Sy(Ω)=∣∣∣∣k=0∑pake−jΩk∣∣∣∣2σw2

R y ( m ) = R y ( n 1 − n 2 ) = E { y ( n 2 ) [ w ( n 1 ) − ∑ k = 1 p a k y ( n 1 − k ) ] } = R w y ( m ) − ∑ k = 1 p a k R y ( m − k ) = σ w 2 h ( − m ) − ∑ k = 1 p a k R y ( m − k ) \begin{aligned} & {{R}_{y}}\left( m \right)={{R}_{y}}\left( {{n}_{1}}-{{n}_{2}} \right)=E\left\{ y\left( {{n}_{2}} \right)\left[ w\left( {{n}_{1}} \right)-\sum\limits_{k=1}^{p}{{{a}_{k}}y\left( {{n}_{1}}-k \right)} \right] \right\} \\ & ={{R}_{wy}}\left( m \right)-\sum\limits_{k=1}^{p}{{{a}_{k}}{{R}_{y}}\left( m-k \right)} \\ & =\sigma _{w}^{2}h\left( -m \right)-\sum\limits_{k=1}^{p}{{{a}_{k}}{{R}_{y}}\left( m-k \right)} \end{aligned} Ry(m)=Ry(n1−n2)=E{y(n2)[w(n1)−k=1∑paky(n1−k)]}=Rwy(m)−k=1∑pakRy(m−k)=σw2h(−m)−k=1∑pakRy(m−k)

同理，将 a 1 = a 2 = ⋯ = a p = 0 {{a}_{1}}={{a}_{2}}=\cdots ={{a}_{p}}=0 a1=a2=⋯=ap=0时，此时为MA（p）模型，分析方法与上述类似，这里不再叙述。

利用MATLAB实现ARMA滤波

MATLAB中有filter函数，其功能是进行一维数字滤波

filter: One-dimensional digital filter. Y = filter(B,A,X) filters the data in vector X with the filter described by vectors A and B to create the filtered data Y. The filter is a “Direct Form II Transposed” implementation of the standard difference equation:

a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)- a(2)*y(n-1) - ... - a(na+1)*y(n-na)If a(1) is not equal to 1, filter normalizes the filter
coefficients by a(1). filter always operates along the first non-singleton dimension,
namely dimension 1 for column vectors and non-trivial matrices,
and dimension 2 for row vectors.[Y,Zf] = filter(B,A,X,Zi) gives access to initial and final
conditions, Zi and Zf, of the delays.  Zi is a vector of length
MAX(LENGTH(A),LENGTH(B))-1, or an array with the leading dimension
of size MAX(LENGTH(A),LENGTH(B))-1 and with remaining dimensions
matching those of X.filter(B,A,X,[],DIM) or filter(B,A,X,Zi,DIM) operates along the
dimension DIM.Tip:  If you have the Signal Processing Toolbox, you can design a
filter, D, using DESIGNFILT.  Then you can use Y = filter(D,X) to
filter your data.

上述是MATLAB中关于filter函数的介绍，可以看出其利用的ARMA模型进行滤波，因此可以将白噪声进行滤波获得色噪声，当然也可以利用其它滤波器进行滤波得到色噪声。
下面给出一个产生色噪声的小例子，样本数为1000。

将白噪声的一些特性画出来可以得到上述的结果。
设计一个低通滤波器（也可以直接指定系数，利用filter函数直接进行滤波），其中滤波器的幅频特性如下：

白噪声经过滤波器后，可以得到色噪声，其结果如下：

可以看出，经过低通滤波后的白噪声的频谱或者是功率谱均不再保持平坦的特性，也就是得到了色噪声。
由于仿真时样本数较少，只能大概看出白噪声的功率谱基本是不变的，如果想得到更好的效果可以增加样本数，这里就不再叙述。

代码如下：

clear;
close all;
clc;
%产生高斯白噪声
N = 1e3;
t=0:(N-1);
w=randn(1,N);
mean_w=mean(w);                      %均值
var_w=var(w);                               %方差
[corr_w,lag]=xcorr(w,'unbiased');   %自相关函数
msv_w = var_w + mean_w.^2;     %均方值
[f1,x_w] = ksdensity(w);               %概率密度
f=0:N-1;
W=fft(w);
fft_w=abs(W);                               %频谱
power_w=W.*conj(W)/1024;          %功率谱密度
figure(1);
subplot(2,2,1);plot(w);
title('高斯白噪声');axis([0 N -5 5]);
subplot(2,2,2);plot(t,mean_w*ones(size(t)));
title('高斯白噪声均值');axis([0 N -2 2]);
subplot(2,2,3);plot(t,var_w*ones(size(t)));
title('高斯白噪声方差');axis([0 N -2 2]);
subplot(2,2,4);plot(t,msv_w*ones(size(t)));
title('高斯白噪声均方值');axis([0 N -2 2]);
figure(2)
subplot(2,2,1);plot(lag,corr_w);
title('高斯白噪声自相关函数');axis([-N N -1 1]);
subplot(2,2,2);plot(x_w,f1);
title('概率密度');
subplot(2,2,3);plot(f,fft_w);
title('高斯白噪声频谱');axis([0 N 0 80]);
subplot(2,2,4);plot(f,power_w);
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;
[N1,Wn]=buttord(2*Wp/W,2*Ws/W,Rp,Rs);
[b,a]=butter(N1,Wn);
[h,f]=freqz(b,a,1000,fs);
figure(3);plot(f,abs(h)); xlabel('f/Hz');ylabel('|H(jf)|');
axis([0 100 0 1.2]);grid on;
title('低通滤波器幅频响应');
%生成高斯色噪声
y=filter(b,a,w);        %滤波产生高斯色噪声
mean_y=mean(y);     %均值
var_y=var(y);               %方差
[corr_y,lag]=xcorr(y,'unbiased');       %自相关函数
msv_y = var_y + mean_y.^2;     %均方值
[f1,x_y]=ksdensity(y);              %概率密度
f=0:N-1;
Y=fft(y);
fft_y=abs(Y);                       %频谱
power_y=Y.*conj(Y)/N;       %功率谱密度
figure(4);
subplot(2,2,1);plot(t,y);
title('高斯色噪声');axis([0 1024 -5 5]);
subplot(2,2,2);plot(t,mean_y*ones(size(t)));
title('高斯色噪声均值');axis([0 N -1 1]);
subplot(2,2,3);plot(t,var_y*ones(size(t)));
title('高斯色噪声方差');axis([0 N -0.5 1.5]);
subplot(2,2,4);plot(t,msv_y*ones(size(t)));
title('高斯色噪声均方值');axis([0 N -0.5 1.5]);
figure(5)
subplot(2,2,1);plot(lag,corr_y);
title('高斯色噪声自相关函数');axis([-N N -0.5 1]);
subplot(2,2,2);plot(x_y,f1);
title('高斯色噪声概率密度');
subplot(2,2,3);plot(f,fft_y);
title('高斯色噪声频谱');axis([0 N 0 100]);
subplot(2,2,4);plot(f,power_y);
title('高斯色噪声功率谱密度');axis([0 N 0 8]);

参考文献：
[1]叶中付. 统计信号处理 [M]. 中国科学技术大学出版社, 2013.