1.白噪声

什么是白噪声

对于时间序列{wt} 若满足下面三个条件，该序列为一个离散的白噪声（white noise）：

每个时间点均值为0：E(wt)=0
每个时间点方差为σ2：Var(wt)=σ2
对于任意k≥1，自相关系ρk=0：Cov(wt,wt+1)=0

使用matlab产生白噪声

fs=10000;                       %采样频率
Ns=10000;                       %采样点数
t=0:1/fs:Ns/fs;
s3=rand(1, Ns+1);%生成一个白噪声数组

图 1‑1白噪声信号

nfft=length(s3);
cxn=xcorr(s3,'unbiased');                  % 计算序列的自相关函数
CXk=fft(cxn,nfft);
Pxx=abs(CXk);
index=0:round(nfft/2-1);
k=index*fs/nfft;
%plot_Pxx=10*log10(Pxx(index+1)); 
plot_Pxx=(Pxx(index+1));

先进行自相关求解，然后进行fft，就能得到功率谱。

图 1‑2白噪声信号功率谱

histogram(s3,25) %直方图

图 1‑3白噪声信号直方图

2.高斯白噪声

什么是高斯噪声

高斯白噪声（white Gaussian noise; WGN）：均匀分布于给定频带上的高斯噪声

如果一个噪声，它的幅度服从高斯分布，而它的功率谱密度又是均匀分布的，则称它为高斯白噪声。
高斯白噪声中的高斯是指：概率分布是正态函数，而白噪声是指：它的二阶矩不相关，一阶矩为常数，是指先后信号在时间上的相关性。这是考察一个信号的两个不同方面的问题。
热噪声和散粒噪声是高斯白噪声。

使用matlab产生白噪声

s3=randn(1, Ns+1);%生成一个高斯白噪声数组

方法与白噪声同理，产生高斯白噪声使用randn函数。

图 2‑1高斯白噪声信号

图 2‑2高斯白噪声信号功率谱

高斯白噪声直方图

图 2‑3高斯白噪声信号直方图

matlab滤波器使用

图 3‑1正弦信号和噪声信号

图 1‑4 红色信号为频率为5Hz的正弦波信号幅值为10，绿色正弦波信号幅值为1，频率为1KHz的信号，黑色信号为白噪声信号。

图 3‑2信号进行相加得到的结果

图 1‑5 为上面三个信号相加得到的结果，1KHz正弦波和白噪声为噪声信号，下面使用FIR滤波器将噪声滤掉。

图 3‑3使用FIR滤波器设计

图 1‑6 使用了matlab中的滤波器设计模块进行设计，设计的采样频率为10K，截至频率为10HZ，阶数为100，得到了上图的结果，最后将它生成为代码，进行程序的编写。

图 3‑4滤波前傅里叶变换

图 3-4可以看出，改信号中包含了10Hz、1KHz以及白噪声的幅值谱。

图 3‑5滤波后信号

图 3-5 可以看出通过滤波器后得到了一个较好的正弦波信号。

图 3‑6滤波后傅里叶变换

图 3-6 可以明显的看出幅值谱已经没有了频率为1K的成分，已经为滤波器滤掉。

clear;
close all;
clc;
f1=10;f2=300;     %原始信号频率
fs=10000;                       %采样频率
Ns=10000;                       %采样点数
t=0:1/fs:Ns/fs;
s1=10*sin(2*pi*f1*t);
s2=sin(2*pi*f2*t);
s3=randn(1, Ns+1);%生成一个白噪声数组
nfft=length(s3);
cxn=xcorr(s3,'unbiased');                  % 计算序列的自相关函数
CXk=fft(cxn,nfft);
Pxx=abs(CXk);
index=0:round(nfft/2-1);
k=index*fs/nfft;
%plot_Pxx=10*log10(Pxx(index+1)); 
plot_Pxx=(Pxx(index+1)); 
figure(1);
subplot(211)
plot(t,s1,'r',t,s2,'g',t,s3,'k');
xlim([0 5000/fs])
title('正弦信号');xlabel('t/s');
y=s1+s2+s3;
set(gcf,'color','w');
subplot(212)
plot(t,y);
title('原始信号');xlabel('t/s');
set(gcf,'color','w');
N=4096;
f=(0:N-1)*fs/N;
yy=abs(fft(y,N)*2/N);
figure(2);
subplot(311)
plot(f(1:N/2),yy(1:N/2));
title('原始信号的FFT');xlabel('f/Hz');
xlim([0 1500]);
h=myfilter;
ylp=filter(h,y);
subplot(312)
plot(t,ylp);
xlim([0.2 0.8]);
title('低通滤波之后的信号')
xlabel('t/s');
yyy=abs(fft(ylp,N)*2/N);
subplot(313)
plot(f(1:N/2),yyy(1:N/2));
xlabel('f/Hz');
title('低通滤波之后的信号频谱');
xlim([0 1500]);
figure(3);
subplot(311)
plot(t,s3);
xlim([0 500/fs])
title('噪声信号');xlabel('t/s');
subplot(312);plot(k,plot_Pxx);title('间接法原始噪声信号的功率谱');     %绘制分贝形式的功率谱
subplot(313);histogram(s3,25);title('噪声直方图');     %直方图

function Hd = myfilter
%MYFILTER 返回离散时间滤波器对象。% MATLAB Code
% Generated by MATLAB(R) 9.8 and DSP System Toolbox 9.10.
% Generated on: 17-Sep-2022 20:01:49% FIR Window Lowpass filter designed using the FIR1 function.% All frequency values are in Hz.
Fs = 10000;  % Sampling FrequencyN             = 100;      % Order
Fc            = 10;       % Cutoff Frequency
flag          = 'scale';  % Sampling Flag
SidelobeAtten = 100;      % Window Parameter% Create the window vector for the design algorithm.
win = chebwin(N+1, SidelobeAtten);% Calculate the coefficients using the FIR1 function.
b  = fir1(N, Fc/(Fs/2), 'low', win, flag);
Hd = dfilt.dffir(b);% [EOF]