Matlab 产生白噪声和有色噪声序列

原文地址 http://blog.sina.com.cn/s/blog_84024a4a01016fmb.html

一、白噪声和有色噪声定义

1.白噪声（white noise）

系统辨识中所用到的数据通常都是含有噪声的。从工程实际出发，这种噪声往往可以视为具有有理谱密度的平稳随机过程。白噪声是一种最简单的随机过程，是有一系列不相关的随机变量组成的理想化随机过程。其自相关函数为dirac函数。

2.有色噪声（colored noise）

理想的白噪声只是一种理论上的抽象，在物理上是很难实现的，现实中并不存在这样的噪声。因而，工程实际中测量数据所包含的噪声往往是有色造势。所谓有色噪声（或相关噪声）是指序列中没一时刻的噪声相关。有色噪声可以看成是由白噪声序列驱动的线性环节的输出。

二、白噪声与有色噪声区别

（1）其实由定义可以看出，白噪声不同时刻是不相关的，自相关函数为脉冲函数；有色噪声则是相关的。

（2）实际测试可以通过测试功率谱来区别，白噪声的功率谱在各频率的值都比较平均，有色噪声则会有较为明显的峰值。

三、具体实例

1.产生有色噪声e(k) = x(k) + 0.5*x(k-1)。其中，x(k)为方差为1的白噪声

clear all; close all;
clc
L=500;  %仿真长度
c = [1 -0.5];
nc = length(c) - 1;
xik=zeros(nc,1);  %白噪声初值
xi=randn(L,1);  %产生均值为0，方差为1的高斯白噪声序列

for k=1:L
    e(k)=c*[xi(k);xik];  %产生有色噪声
    %数据更新
    for i=nc:-1:2
        xik(i)=xik(i-1);
    end
    xik(1)=xi(k);
end

subplot(2,1,1);
plot(xi);
xlabel('k');ylabel('噪声幅值');title('白噪声序列');
subplot(2,1,2);
plot(e);
xlabel('k');ylabel('噪声幅值');title('有色噪声序列');

%测试功率谱

[y1,f1] = Spectrum_Calc(xi',512);
p1 = 1/L * y1.*conj(y1);

figure(2)
subplot(211)
plot(f1,p1)

[y2,f2] = Spectrum_Calc(e,512);
p2 = 1/L * y2.*conj(y2);
subplot(212)
plot(f2,p2)

 

运行结果：

2：

 

clear all; close all;
L=500;  %仿真长度
d=[1 -1.5 0.7 0.1]; c=[1 0.5 0.2];  % 分子分母多项式系数
nd=length(d)-1 ;nc=length(c)-1;   %阶次
xik=zeros(nc,1);  %白噪声初值
ek=zeros(nd,1);
xi=randn(L,1);  %产生均值为0，方差为1的高斯白噪声序列

for k=1:L
    e(k)=-d(2:nd+1)*ek+c*[xi(k);xik];  %产生有色噪声
    %数据更新
    for i=nd:-1:2
        ek(i)=ek(i-1);
    end
    ek(1)=e(k);
    for i=nc:-1:2
        xik(i)=xik(i-1);
    end
    xik(1)=xi(k);
end
subplot(2,1,1);
plot(xi);
xlabel('k');ylabel('噪声幅值');title('白噪声序列');
subplot(2,1,2);
plot(e);
xlabel('k');ylabel('噪声幅值');title('有色噪声序列');

%测试功率谱

[y1,f1] = Spectrum_Calc(xi',512);
p1 = 1/L * y1.*conj(y1);

figure(2)
subplot(211)
plot(f1,p1)

[y2,f2] = Spectrum_Calc(e,512);
p2 = 1/L * y2.*conj(y2);
subplot(212)
plot(f2,p2)

这个功率谱比例1效果要好很多。

四、程序及结果说明：

1.Spectrum_Calc函数为快速傅里叶变换函数，具体见博文

http://blog.sina.com.cn/s/blog_84024a4a01015rez.html

2.例2中提供了一个线性系统，将系统转换为差分方程形式可能更有助于程序的理解。系统的等价差分方程为：

e(k) – 1.5e(k-1) + 0.7e(k-2) + 0.1e(k-3) = xi(k) + 0.5xi(k-1) + 0.2xi(k-2)