17、matlab实现均值滤波、中值滤波、Butterworth滤波和线性相位FIR滤波

1、创建信号

1）创建正余弦信号、噪声信号和混合信号

原始正余弦信号公式：Signal1= sin(2*pi*20* t) + sin(2*pi*40* t) + sin(2*pi*60* t)

高斯分布的白噪声：NoiseGauss= [randn(1,2000)]

均匀分布的白噪声：[rand(1,2000)]

正余弦信号、噪声信号和混合信号代码：

N = 2000; %采样点数
Fs = 2000; %采样频率
t = 0:1 / Fs:1 - 1 / Fs; %时间序列
Signal1= sin(2*pi*20* t) + sin(2*pi*40* t) + sin(2*pi*60* t);
NoiseGauss= [randn(1,2000)]; %高斯分部白噪声
NoiseWhite= [rand(1,2000)]; %后100点均匀分布白噪声
NoiseSignal1= Signal1+NoiseGauss; %设计混合信号1
NoiseSignal2= Signal1+NoiseWhite; %设计混合信号2
subplot(5,1,1);
plot(Signal1);
title('正余弦信号');
subplot(5,1,2);
plot(NoiseGauss);
title('高斯噪声');
subplot(5,1,3);
plot(NoiseWhite);
title('白噪声');
subplot(5,1,4);
plot(NoiseSignal1);
title('混合高斯噪声信号');
subplot(5,1,5);
plot(NoiseSignal2);
title('混合白噪声信号');

试图效果：

2）创建方波信号、噪声及混合信号

原始方波信号公式：

Signal2=[2*ones(1,50),zeros(1,50),-1*ones(1,100),zeros(1,50),-2*ones(1,50),zeros(1,50),1*ones(1,100),zeros(1,50),2*ones(1,50),zeros(1,50),-1*ones(1,100),zeros(1,50),-2*ones(1,50),zeros(1,50),1*ones(1,100),zeros(1,50),2*ones(1,50),zeros(1,50),-1*ones(1,100),zeros(1,50),-2*ones(1,50),zeros(1,50),1*ones(1,100),zeros(1,50),2*ones(1,50),zeros(1,50),-1*ones(1,100),zeros(1,50),-2*ones(1,50),zeros(1,50),1*ones(1,100),zeros(1,50),]

高斯分布的白噪声：NoiseGauss= [randn(1,2000)]

均匀分布的白噪声：[rand(1,2000)]

方波信号、噪声及混合信号代码

N = 2000; %采样点数
Fs = 2000; %采样频率
t = 0:1 / Fs:1 - 1 / Fs; %时间序列
NoiseGauss= [randn(1,2000)]; %高斯分部白噪声
NoiseWhite= [rand(1,2000)]; %后100点均匀分布白噪声
Signal2=[2*ones(1,50),zeros(1,50),-1*ones(1,100),zeros(1,50),-2*ones(1,50),zeros(1,50),1*ones(1,100),zeros(1,50),2*ones(1,50),zeros(1,50),-1*ones(1,100),zeros(1,50),-2*ones(1,50),zeros(1,50),1*ones(1,100),zeros(1,50),2*ones(1,50),zeros(1,50),-1*ones(1,100),zeros(1,50),-2*ones(1,50),zeros(1,50),1*ones(1,100),zeros(1,50),2*ones(1,50),zeros(1,50),-1*ones(1,100),zeros(1,50),-2*ones(1,50),zeros(1,50),1*ones(1,100),zeros(1,50),];
NSignal1= Signal2+NoiseGauss; %设计混合信号1
NSignal2= Signal2+NoiseWhite; %设计混合信号2
subplot(5,1,1);
plot(Signal2);
title('方波信号');
subplot(5,1,2);
plot(NoiseGauss);
title('高斯噪声');
subplot(5,1,3);
plot(NoiseWhite);
title('白噪声');
subplot(5,1,4);
plot(NSignal1);
title('方波混合高斯噪声信号');
subplot(5,1,5);
plot(NSignal2);
title('方波混合白噪声信号');

试图效果：

2、均值滤波filter()函数

语法：y = filter(b,a,x) 使用由分子和分母系数 b 和 a 定义的有理传递函数对输入数据 x 进行滤波。

1）正余弦混合噪声信号均值滤波

代码：

N = 2000; %采样点数
Fs = 2000; %采样频率
t = 0:1 / Fs:1 - 1 / Fs; %时间序列
Signal1= sin(2*pi*20* t) + sin(2*pi*40* t) + sin(2*pi*60* t);
NoiseGauss= [randn(1,2000)]; %高斯分部白噪声
NoiseWhite= [rand(1,2000)]; %后100点均匀分布白噪声
NoiseSignal1= Signal1+NoiseGauss; %设计混合信号1
NoiseSignal2= Signal1+NoiseWhite; %设计混合信号2
b=[1 1 1 1 1 1]/6;
Signal_Filter1 = filter(b,1,NoiseSignal1);
Signal_Filter2 = filter(b,1,NoiseSignal2);
figure(1);
subplot(3,2,1);               
plot(Signal1);
title('原始信号1');
subplot(3,2,3);               
plot(NoiseSignal1);
title('高斯混合信号1');
subplot(3,2,5);  
plot(Signal_Filter1);
title('均值滤波后高斯混合信号1');
subplot(3,2,2);               
plot(Signal1);
title('原始信号1');
subplot(3,2,4);               
plot(NoiseSignal2);
title('白噪声混合信号1');
subplot(3,2,6);  
plot(Signal_Filter2);
title('均值滤波后白噪声混合信号1');

视图效果：竖直方向三幅图进行对比

2）方波混合噪声信号均值滤波

代码

N = 2000; %采样点数
Fs = 2000; %采样频率
t = 0:1 / Fs:1 - 1 / Fs; %时间序列
NoiseGauss= [randn(1,2000)]; %高斯分部白噪声
NoiseWhite= [rand(1,2000)]; %后100点均匀分布白噪声
Signal2=[2*ones(1,50),zeros(1,50),-1*ones(1,100),zeros(1,50),-2*ones(1,50),zeros(1,50),1*ones(1,100),zeros(1,50),2*ones(1,50),zeros(1,50),-1*ones(1,100),zeros(1,50),-2*ones(1,50),zeros(1,50),1*ones(1,100),zeros(1,50),2*ones(1,50),zeros(1,50),-1*ones(1,100),zeros(1,50),-2*ones(1,50),zeros(1,50),1*ones(1,100),zeros(1,50),2*ones(1,50),zeros(1,50),-1*ones(1,100),zeros(1,50),-2*ones(1,50),zeros(1,50),1*ones(1,100),zeros(1,50),];
NSignal1= Signal2+NoiseGauss; %设计混合信号1
NSignal2= Signal2+NoiseWhite; %设计混合信号2
b=[1 1 1 1 1 1]/6;
Signal_Filter3 = filter(b,1,NSignal1);
Signal_Filter4 = filter(b,1,NSignal2);
figure(2);
subplot(3,2,1);               
plot(Signal2);
title('原始信号2');
subplot(3,2,3);               
plot(NSignal1);
title('高斯混合信号2');
subplot(3,2,5);  
plot(Signal_Filter3);
title('均值滤波后高斯混合信号2');
subplot(3,2,2);               
plot(Signal2);
title('原始信号2');
subplot(3,2,4);               
plot(NSignal2);
title('白噪声混合信号2');
subplot(3,2,6);  
plot(Signal_Filter4);
title('均值滤波后白噪声混合信号2');

视图效果：

3、中值滤波：medfilt1()函数

语法：medfilt1(x,n)对x应用一个n阶一维中值滤波器

1）正余弦混合噪声信号中值滤波

代码

N = 2000; %采样点数
Fs = 2000; %采样频率
t = 0:1 / Fs:1 - 1 / Fs; %时间序列
Signal1= sin(2*pi*20* t) + sin(2*pi*40* t) + sin(2*pi*60* t);
NoiseGauss= [randn(1,2000)]; %高斯分部白噪声
NoiseWhite= [rand(1,2000)]; %后100点均匀分布白噪声
NoiseSignal1= Signal1+NoiseGauss; %设计混合信号1
NoiseSignal2= Signal1+NoiseWhite; %设计混合信号2
Signal_Filter1=medfilt1(NoiseSignal1,11);
Signal_Filter2=medfilt1(NoiseSignal2,11);
figure(1);
subplot(3,2,1);               
plot(Signal1);
title('原始信号1');
subplot(3,2,3);               
plot(NoiseSignal1);
title('高斯混合信号1');
subplot(3,2,5);  
plot(Signal_Filter1);
title('中值滤波后高斯混合信号1');
subplot(3,2,2);               
plot(Signal1);
title('原始信号1');
subplot(3,2,4);               
plot(NoiseSignal2);
title('白噪声混合信号1');
subplot(3,2,6);  
plot(Signal_Filter2);
title('中值滤波后白噪声混合信号1');

视图效果

2）方波混合噪声信号中值滤波

代码

N = 2000; %采样点数
Fs = 2000; %采样频率
t = 0:1 / Fs:1 - 1 / Fs; %时间序列
NoiseGauss= [randn(1,2000)]; %高斯分部白噪声
NoiseWhite= [rand(1,2000)]; %后100点均匀分布白噪声
Signal2=[2*ones(1,50),zeros(1,50),-1*ones(1,100),zeros(1,50),-2*ones(1,50),zeros(1,50),1*ones(1,100),zeros(1,50),2*ones(1,50),zeros(1,50),-1*ones(1,100),zeros(1,50),-2*ones(1,50),zeros(1,50),1*ones(1,100),zeros(1,50),2*ones(1,50),zeros(1,50),-1*ones(1,100),zeros(1,50),-2*ones(1,50),zeros(1,50),1*ones(1,100),zeros(1,50),2*ones(1,50),zeros(1,50),-1*ones(1,100),zeros(1,50),-2*ones(1,50),zeros(1,50),1*ones(1,100),zeros(1,50),];
NSignal1= Signal2+NoiseGauss; %设计混合信号1
NSignal2= Signal2+NoiseWhite; %设计混合信号2
Signal_Filter3=medfilt1(NSignal1,11);
Signal_Filter4=medfilt1(NSignal2,11);
figure(2);
subplot(3,2,1);               
plot(Signal2);
title('原始信号2');
subplot(3,2,3);               
plot(NSignal1);
title('高斯混合信号2');
subplot(3,2,5);  
plot(Signal_Filter3);
title('中值滤波后高斯混合信号2');
subplot(3,2,2);               
plot(Signal2);
title('原始信号2');
subplot(3,2,4);               
plot(NSignal2);
title('白噪声混合信号2');
subplot(3,2,6);  
plot(Signal_Filter4);
title('中值滤波后白噪声混合信号2');

视图效果

4、 Butterworth低通滤波：巴特沃斯滤波butter()函数

语法：[b,a] = butter(n,Wn) 返回归一化截止频率为 Wn 的 n 阶低通数字巴特沃斯滤波器的传递函数系数。

1）正余弦混合噪声信号Butterworth低通滤波
代码

N = 2000; %采样点数
Fs = 2000; %采样频率
t = 0:1 / Fs:1 - 1 / Fs; %时间序列
Signal1= sin(2*pi*20* t) + sin(2*pi*40* t) + sin(2*pi*60* t);
NoiseGauss= [randn(1,2000)]; %高斯分部白噪声
NoiseWhite= [rand(1,2000)]; %后100点均匀分布白噪声
NoiseSignal1= Signal1+NoiseGauss; %设计混合信号1
NoiseSignal2= Signal1+NoiseWhite; %设计混合信号2
Wc=2*250/Fs;%截止频率 50Hz
[b,a]=butter(4,Wc);%Butterworth滤波
Signal_Filter1=filter(b,a,NoiseSignal1);
Signal_Filter2=filter(b,a,NoiseSignal2);
figure(1);
subplot(3,2,1);               
plot(Signal1);
title('原始信号1');
subplot(3,2,3);               
plot(NoiseSignal1);
title('高斯混合信号1');
subplot(3,2,5);  
plot(Signal_Filter1);
title('巴特沃斯低通滤波后高斯混合信号1');
subplot(3,2,2);               
plot(Signal1);
title('原始信号1');
subplot(3,2,4);               
plot(NoiseSignal2);
title('白噪声混合信号1');
subplot(3,2,6);  
plot(Signal_Filter2);
title('巴特沃斯低通滤波后白噪声混合信号1');

视图效果

2）方波混合噪声信号Butterworth低通滤波

代码

N = 2000; %采样点数
Fs = 2000; %采样频率
t = 0:1 / Fs:1 - 1 / Fs; %时间序列
NoiseGauss= [randn(1,2000)]; %高斯分部白噪声
NoiseWhite= [rand(1,2000)]; %后100点均匀分布白噪声
Signal2=[2*ones(1,50),zeros(1,50),-1*ones(1,100),zeros(1,50),-2*ones(1,50),zeros(1,50),1*ones(1,100),zeros(1,50),2*ones(1,50),zeros(1,50),-1*ones(1,100),zeros(1,50),-2*ones(1,50),zeros(1,50),1*ones(1,100),zeros(1,50),2*ones(1,50),zeros(1,50),-1*ones(1,100),zeros(1,50),-2*ones(1,50),zeros(1,50),1*ones(1,100),zeros(1,50),2*ones(1,50),zeros(1,50),-1*ones(1,100),zeros(1,50),-2*ones(1,50),zeros(1,50),1*ones(1,100),zeros(1,50),];
NSignal1= Signal2+NoiseGauss; %设计混合信号1
NSignal2= Signal2+NoiseWhite; %设计混合信号2
Wc=2*250/Fs;%截止频率 50Hz
[b,a]=butter(4,Wc);%Butterworth滤波
Signal_Filter3=filter(b,a,NSignal1);
Signal_Filter4=filter(b,a,NSignal2);
figure(2);
subplot(3,2,1);               
plot(Signal2);
title('原始信号2');
subplot(3,2,3);               
plot(NSignal1);
title('高斯混合信号2');
subplot(3,2,5);  
plot(Signal_Filter3);
title('巴特沃斯低通滤波后高斯混合信号2');
subplot(3,2,2);               
plot(Signal2);
title('原始信号2');
subplot(3,2,4);               
plot(NSignal2);
title('白噪声混合信号2');
subplot(3,2,6);  
plot(Signal_Filter4);
title('巴特沃斯低通滤波后白噪声混合信号2');

视图效果

4、线性相位FIR滤波：firls()函数

语法：firls(n,f,a)：返回包含n阶FIR滤波器的n+1个系数的行向量b。所得滤波器的频率和幅度特性与向量f和a给出的特性相匹配。

1）正余弦混合噪声信号线性相位FIR滤波

代码

N = 2000; %采样点数
Fs = 2000; %采样频率
t = 0:1 / Fs:1 - 1 / Fs; %时间序列
Signal1= sin(2*pi*20* t) + sin(2*pi*40* t) + sin(2*pi*60* t);
NoiseGauss= [randn(1,2000)]; %高斯分部白噪声
NoiseWhite= [rand(1,2000)]; %后100点均匀分布白噪声
NoiseSignal1= Signal1+NoiseGauss; %设计混合信号1
NoiseSignal2= Signal1+NoiseWhite; %设计混合信号2
F  =  0:0.05:0.95;%频率点对的向量
A  =  [1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];%与F的参数对应
b  =  firls(2,F,A);%20为阶次
Signal_Filter1 = filter(b,1,NoiseSignal1);
Signal_Filter2 = filter(b,1,NoiseSignal2);
figure(1);
subplot(3,2,1);               
plot(Signal1);
title('原始信号1');
subplot(3,2,3);               
plot(NoiseSignal1);
title('高斯混合信号1');
subplot(3,2,5);  
plot(Signal_Filter1);
title('FIR滤波后高斯混合信号1');
subplot(3,2,2);               
plot(Signal1);
title('原始信号1');
subplot(3,2,4);               
plot(NoiseSignal2);
title('白噪声混合信号1');
subplot(3,2,6);  
plot(Signal_Filter2);
title('FIR滤波后白噪声混合信号1');

视图效果

2）方波混合噪声信号线性相位FIR滤波

代码

N = 2000; %采样点数
Fs = 2000; %采样频率
t = 0:1 / Fs:1 - 1 / Fs; %时间序列
NoiseGauss= [randn(1,2000)]; %高斯分部白噪声
NoiseWhite= [rand(1,2000)]; %后100点均匀分布白噪声
Signal2=[2*ones(1,50),zeros(1,50),-1*ones(1,100),zeros(1,50),-2*ones(1,50),zeros(1,50),1*ones(1,100),zeros(1,50),2*ones(1,50),zeros(1,50),-1*ones(1,100),zeros(1,50),-2*ones(1,50),zeros(1,50),1*ones(1,100),zeros(1,50),2*ones(1,50),zeros(1,50),-1*ones(1,100),zeros(1,50),-2*ones(1,50),zeros(1,50),1*ones(1,100),zeros(1,50),2*ones(1,50),zeros(1,50),-1*ones(1,100),zeros(1,50),-2*ones(1,50),zeros(1,50),1*ones(1,100),zeros(1,50),];
NSignal1= Signal2+NoiseGauss; %设计混合信号1
NSignal2= Signal2+NoiseWhite; %设计混合信号2
F  =  0:0.05:0.95;%频率点对的向量
A  =  [1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];%与F的参数对应
b  =  firls(2,F,A);%20为阶次
Signal_Filter3 = filter(b,1,NSignal1);
Signal_Filter4 = filter(b,1,NSignal2);
figure(2);
subplot(3,2,1);               
plot(Signal2);
title('原始信号2');
subplot(3,2,3);               
plot(NSignal1);
title('高斯混合信号2');
subplot(3,2,5);  
plot(Signal_Filter3);
title('FIR滤波后高斯混合信号2');
subplot(3,2,2);               
plot(Signal2);
title('原始信号2');
subplot(3,2,4);               
plot(NSignal2);
title('白噪声混合信号2');
subplot(3,2,6);  
plot(Signal_Filter4);
title('FIR滤波后白噪声混合信号2');