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);%生成一个白噪声数组
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‑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]