本文介绍了MATLAB中数字滤波技术的应用，特别是中值滤波和算术平均值滤波。中值滤波在去除椒盐噪声方面效果显著，适用于温度、液位等变化缓慢的参数；算术平均值滤波则适用于平滑处理信号。文中提供了一个原创的中值滤波函数Median_Filter，可用于二值图像、RGB图像或函数图像，并展示了其在不同场景下的应用。同时，还介绍了算术平均值滤波函数Mean_Filter，用于对比中值滤波的效果。
阅读前请注意：
1. MATLAB深入学习是在MATLAB学习笔记的基础上的深度学习，学习这一部分内容需要有较夯实的MATLAB基础，不建议在未掌握基础知识的情况下学习本内容。博客内容由 @K2SO4钾 撰写、编辑，发布于 @K2SO4钾 的个人投稿与华中师范大学HelloWorld程序设计协会CSDN官方账号 @CCNU_HelloWorld 。注意，如需转载需得到作者 @K2SO4钾 的同意与授权！
2. 内容中的代码、示例均为作者原创，内容中的引用、参考等非原创内容均会以参考文献的方式在文末标出。
3. 笔记中示例用MATLAB版本为MATLAB R2019a。
4. 请谅解笔记可能会出现的错误，欢迎指正讨论；由于MATLAB更新导致的旧代码无法使用的情况，也欢迎讨论交流。
文章目录
S1.0 MATLAB学习笔记：传送门汇总
S1.0.1 MATLAB学习笔记：传送门汇总
S1.0.2 MATLAB拓展学习：传送门汇总
MATLAB拓展学习T2：程序性能优化技术
MATLAB拓展学习T3：histogram函数详解
MATLAB拓展学习T4：数据导入与查表技术
S1.0.3 MATLAB深入学习：传送门汇总
MATLAB深入学习S1：简单数字滤波技术的MATLAB实现（1）
MATLAB深入学习S2：无限脉冲响应数字滤波器
MATLAB深入学习S3：简单数字滤波技术的MATLAB实现
MATLAB深入学习S4：元胞自动机原理及MATLAB实现
MATLAB深入学习S5：图像处理技术
S1.0.4 MATLAB应用实例：传送门汇总
MATLAB应用实例Y1：Floyd算法
S1.1 数字滤波概述
对于一个信号，我们常常使用滤波技术保留对我们有益的波段，去除无益的波段、噪声和干扰。但实际在工业领域中，往往会有较多的强干扰源。在环境如温度、电场、磁场、湿度等作用下，加之设备本身造成的干扰，往往会产生大量随机出现的干扰信号。干扰信号不仅会影响设备对数据的收集，在工业过程控制系统中也会影响受控对象和整个系统工作的可靠性。
简单句几个例子：体重秤有时会因为称重者在秤上的重心变化而剧烈跳动，电脑手机屏幕在附近有强磁场时会出现雪花点，录音设备在录音时常常会收录一些刺耳的杂音，接收某一信号时同时也会接收附近干扰源的干扰信号等。
数字滤波技术，可以使微型计算机系统在采集模拟信号的同时，对收集到的序列进行加工处理（平滑加工等），使得有用信号占收集到信号的大部分比例，尽可能减少环境、设备本身等因素造成的干扰和噪声，保障系统工作的可靠性。相较于传统的模拟滤波，数字滤波可以在微型计算机中设定一段处理程序，以使得该计算机系统在收集信号的同时对信号进行处理。
数字滤波处理同样在计算机图像处理方面起重要作用。数字滤波技术可以使人们用程序的方式加工图片或视频（视频本身也是由多帧图片组成）使得其符合需求。很多情况下，数字滤波技术、图像处理技术能让人们更好地认识世界、处理问题，例如图像增强技术、图像去噪技术能让图片的关键信息更加清晰可见，图像分类识别技术能更好地为人工智能领域服务。数字滤波技术与现代科技发展密不可分。
数字滤波技术与模拟滤波技术相比有以下优点：1
- 数字滤波可以以程序的方式放置在数据处理和控制算法之间，减少硬件设备的消耗；
- 尽可能减少了硬件设备对系统带来的不稳定性（硬件设备本身会有一定损耗）；
- 模拟滤波器各通道通常是专用的，数字滤波器则通常可多通道共享，降低成本；
- 模拟滤波器受电容容量限制，不可能像数字滤波器一样能滤波很低频率的信号；
- 使用更灵活便捷，可进行多次滤波、多种不同滤波方式、多样化且便于调节的滤波参数。
S1.2 常用的简单数字滤波技术
S1.2.1 中值滤波
S1.2.1.1 中值滤波简介
中值滤波对于去掉偶然因素引起的波动或采样器不稳定造成的误差所引起的误差干扰比较有效，对温度、液位等变化缓慢的被控参数有良好的滤波效果，但是对流量、速度等变化较快的参数不宜使用。中值滤波法对过滤椒盐噪声效果显著。2
这是本人在打2021年数模国赛的时候，从华师数学与统计学院6号楼朝向物院9号楼拍摄的照片。理科大楼上有蓝天白云，下有茂林掩映，十分好看。
t = imread('D:\Huawei Share\OneHop\IMG_20210912_122924.jpg');
用imread函数将这张图片导入MATLAB，会在工作区中看到一个 3000 ∗ 4000 ∗ 3 3000*4000*3 3000∗4000∗3的uint类型的矩阵。这个矩阵即描述了这张照片所有像素点的信息。 3000 ∗ 4000 3000*4000 3000∗4000表示这是一个长、宽分别为3000像素点、4000像素点的图片，而三维矩阵中的第三维度，每一个维度依次存储R、G、B值，描述了每一个像素点具体的颜色。矩阵中每一个元素值均为0~255的整数。
如果将这个RGB图像转换为灰度图像，就会在工作区中看到一个 3000 ∗ 4000 3000*4000 3000∗4000的uint类型的矩阵。这个矩阵中每一个元素描述了这张灰度照片每一个像素点灰度值，灰度值越大越接近于白色，灰度值越小越接近于黑色。
t = rgb2gray(t);
现在，我们向RGB图片里添加一些椒盐噪声：
t = imnoise(t,'salt & pepper',0.1); % 加入密度为0.1的椒盐噪声
所谓的 “椒盐噪声”，就是RGB图像中出现的一些随机颜色的杂点，也是灰度图像中出现的一些随机灰度值的杂点，也称“脉冲噪声”，是图像处理中常见的噪声。下图分别是有椒盐噪声的RGB图像、有椒盐噪声的灰度图像在放大之后的杂点，你看下图灰度图像中的椒盐噪声，像不像你在案板上随机撒下的盐和黑胡椒？
椒盐噪声的成因可能是影像讯号受到突如其来的强烈干扰而产生、类比数位转换器或位元传输错误等。例如失效的感应器导致像素值为最小值，饱和的感应器导致像素值为最大值。常用的去除这种噪声的有效手段是使用中值滤波器。3
一个函数，就比如说是 s i n sin sin函数，如果在 s i n sin sin函数所有采样点的随机处添加一些随机的扰动，如何尽可能消除这些扰动呢，中值滤波法也可以纳入考虑。
clc,clear all,close all
x = 0:0.01:2*pi;
y = sin(x);
subplot(121)
plot(x,y,'b-','LineWidth',1)
title('sin函数')
noise = randi([1,length(x)],1,100);
y(1,noise) = y(1,noise)+randi([-500,500],1,length(noise))/1000;
subplot(122)
plot(x,y,'b-','LineWidth',1)
title('添加随机扰动')
S1.2.1.2 中值滤波的基本原理
MATLAB自带有中值滤波函数medfilt1、medfilt2，分别对一维序列、二维数据进行中值滤波。详细使用方法见帮助文档，这里主要讲解中值滤波的基本原理和函数创建。
中值滤波的思想很简单。对于函数的某一点，取该点附近的几个点作为“取样窗口”，将该点的值用取样窗口值的中位数代替。对于图像中的某一个像素，取该像素点附近的几个像素点作为“取样窗口”，将该像素点的灰度值或是RGB值用取样窗口的灰度值或RGB值的中位数代替。 但是实际应用上，常常面临以下问题：
1. 取样窗口的大小如何选取？
取样窗口的大小选取不同，对于一些函数或图像，最后滤波的结果差异也较大。
窗口大小的选取常考虑的因素有：噪声的密度、序列或图像的大小、噪声本身的长度和宽度（如果窗口长宽小于噪声长宽，那么中值滤波很难除去噪声）、噪声变化的速率、去噪的成本和算法的运行速度。通常，一个序列经过若干次中值滤波后（窗口任意大，但是每一次中值滤波的窗口大小应相同），最终会得到一个稳定的序列。这个序列被称为 “根序列”。根序列对中值滤波窗口大小的选取有重要意义，这里不展开。
通常，人们会根据经验还有若干次实验结果来选定最佳的取样窗口大小。
2. 取样窗口的形状如何选取？
同样的道理，人们会根据经验还有若干次实验结果来选定最佳的取样窗口形状。常见的窗口形状有：正方形，矩形，十字形。下面给出的示例函数也会采用这三种图形。
3. 边缘如何处理？
对于任何的图像处理技术，边缘都是很难处理。由于中值滤波算法需要选取某一个像素点周围的像素点，取中位数代替原像素点（如果是一个函数序列的话就选取某一点附近的点，取中位数代替原点），在边缘的像素点（或点）就无法通过这种方法来处理。
边缘的处理方法常用的有：Canny算子边缘检测法、Laplace算子边缘检测法等，这里不进行展开，下面的中值滤波函数暂时不会给出边缘检测的部分。关于边缘检测，可参考这一篇博客：图像滤波（中值、均值）及边缘检测（matlab实现）。
除这些之外，实际应用中的中值滤波法肯定还需要考虑更多问题。在MATLAB中实现中值滤波会相较于用汇编来实现更加便捷，但是也应该注意下面一些要点：
MATLAB中图像对应的矩阵是uint8类型，而在对数据进行处理时，数据的类型常常会变为double。最终应通过强制转化将矩阵转化为uint8类型；
对于一个任意形状的取样窗口，都应有一个中心点。窗口内所有元素的中位数将代替中心点的值。窗口的中心点到窗口行边缘的距离为“横臂”，即代码中c_arm；窗口的中心点到窗口列边缘的距离为“纵臂”，即代码中r_arm；
由于窗口本身的大小，图像或序列的边界无法处理（至少在讨论这一小节时我们暂时不对边缘进行处理），因此用原图的边缘值或原序列的边缘值取代即可。此外还需要注意，在循环各点进行中值处理的时候，不要把边缘值也处理了（否则窗口就会延伸到图像外，会报错）。
下面给出了一个原创的、较为完整的中值滤波函数Median_Filter，其详细的功能及使用方法请查看这个函数的帮助文档。
function [Output1, Output2] = Median_Filter(Method, Input1, Input2, Window, Window_type, Edge_Check, Drawing)
%
% Return the median filtering result for binary image, RGB image
% or function image in rectangular coordinate system.
%
% Median_Filter函数用于返回对于二值图像、RGB图像或直角坐标系下函数图像的中值滤波结果
%
% 输入参数Method：取值可为1、2、3，其它值的输入会报错
% > Method = 1 时对直角坐标系下的函数图像进行中值滤波
% > Method = 2 时对二值图像进行中值滤波
% > Method = 3 时对RGB图像进行中值滤波
% 输入参数Input：
% > 若Method = 1，则Input1、Input2分别输入函数对应的横纵坐标取值（均为行向量），
% 或者Input1输入空矩阵，Input2输入纵坐标取值，该函数将会从1开始自动填充横坐标
% > 若Method = 2，则Input1输入二值图像的二维矩阵，Input2输入空矩阵
% > 若Method = 3，则Input1输入RGB图像的三维矩阵
% 输入参数Window：中值滤波的窗口大小，（其矩阵中的元素应均为正整数，若输入为双精度数则会向下取整）
% 输入参数Window_type：窗口类型
% > 若Window_type = 'default'（即窗口为一个行向量）
% ·Method = 1时，Window取一个小于Input2向量长度的整数，表示窗口长度（应为大于1的奇数）
% ·Method取其它值时会报错
% > 若Window_type = 'square'（即窗口为正方形）
% ·Method = 1时会报错
% ·若Method = 2或3，Window取一个小于Input1矩阵行宽的整数，表示窗口边长（应为大于1的奇数）
% > 若Window_type = 'rectangle'（即窗口为矩形）
% ·Method = 1时会报错
% ·若Method = 2或3，Window取一个小于Input1矩阵行宽的1*2向量，分别表示窗口长、宽（均应为大于1的奇数）
% > 若Window_type = 'cross'（即窗口为十字形）
% ·Method = 1时会报错
% ·若Method = 2或3，Window取一个小于Input1矩阵行宽的1*2向量，分别表示十字窗口长、宽（均应为大于1的奇数）
% > 若Method = 3，如果需要对R、G、B进行不同窗口大小的中值滤波，请输入3*2或3*1的矩阵（视Window_type而定），
% 每一行分别表示对R、G、B进行中值滤波处理的窗口大小
% 输入参数Edge_Check：边缘检测方式
% > 'none'（默认）：不进行边缘检测
% 输入参数Drawing：是否在函数内进行绘图或展示图像
% > 'off'（默认）：不展示图像
% > 'on'：展示图像
%
% 输出参数Output1：Method = 1时，输出处理结果的横坐标向量；Method = 2或3时输出图像对应的矩阵
% 输出参数Output2：Method = 1时，输出处理结果的纵坐标向量；Method = 2或3时Output2为空
%
% 常用调用方式：
% [Output1, Output2] = Median_Filter(Method, Input1, Input2, Window, Window_type)
% [Output1, Output2] = Median_Filter(Method, Input1, Input2, Window, Window_type, Edge_Check, Drawing)
%
% 编写于2021年10月6日20:32 程序员：K2SO4钾（华中师范大学）
%% 初判断错误与警告
% 设定初始参数
if isempty(Input1)
r_length = size(Input2,1);
c_length = size(Input2,2);
else
r_length = size(Input1,1);
c_length = size(Input1,2);
end
if nargin == 4
Window_type = 'default';
Edge_Check = 'none';
Drawing = 'off';
elseif nargin == 5
Edge_Check = 'none';
Drawing = 'off';
elseif nargin == 6
if strcmp(Edge_Check,'on') || strcmp(Edge_Check,'off')
Drawing = Edge_Check;
Edge_Check = 'none';
else
Drawing = 'off';
end
end
% 输入输出参数规范
if nargin >= 8
error('输入参数数量太多')
elseif nargin <= 3
error('输入参数数量太少')
end
if nargout >= 3
error('输出参数数量太多')
end
% Method参数规范
Method = floor(Method);
if Method >= 4 || Method < 1
error('Method的值规定为1、2、3，其它向下取整后在规定值之外的输入值将会报错')
end
% Method = 1时规定Input1、Input2应为等长的行向量
if Method == 1
if ~isempty(Input1)
if size(Input1,1) ~= 1 || size(Input2,1) ~= 1
error('Input1不为空时，Input1、Input2应为等长的行向量')
end
if length(Input1) ~= length(Input2)
error('Input1不为空时，Input1、Input2应为等长的行向量')
end
else
if size(Input2,1) ~= 1
error('Input2应为行向量')
end
end
end
% Method = 2、3时规定Input1的长宽必须在3*3以上
if Method == 2 || Method == 3
if size(Input1,1) <= 3 || size(Input1,2) <= 3
error('Method = 2、3时规定Input1的长宽必须在3*3以上')
end
if ~isempty(Input2)
warning('将忽略Input2的输入')
end
end
% 窗口大小应设定为正整数、大于1的奇数，且窗口长度或长宽应小于等于函数因变量长度或图像长宽
Window = floor(Window);
if sum(sum(Window < 3)) ~= 0
error('窗口大小应设定为正整数、大于1的奇数')
end
if Method == 1
if Window > c_length
error(error('窗口长度或长宽应小于等于函数因变量长度或图像长宽'))
end
elseif strcmp(Window_type,'square')
if max(Window(:,1)) > r_length || max(Window(:,1)) > c_length
error(error('窗口长度或长宽应小于等于函数因变量长度或图像长宽'))
end
else
if max(Window(:,1)) > r_length || max(Window(:,2)) > c_length
error('窗口长度或长宽应小于等于函数因变量长度或图像长宽')
end
end
% Method与Window_type搭配不符合规定。详情请键入 help Median_Filter
if Method == 1
if ~(strcmp(Window_type,'default'))
error('Method与Window_type搭配不符合规定。详情请键入 help Median_Filter')
end
else
if ~((strcmp(Window_type,'square')) || (strcmp(Window_type,'rectangle')) || (strcmp(Window_type,'cross')))
error('Method与Window_type搭配不符合规定。详情请键入 help Median_Filter')
end
end
% Window矩阵的大小不符合规定。按照规则输入参数Window，详情请键入 help Median_Filter
if Method == 1
if size(Window,1)*size(Window,2) ~= 1
error('Window矩阵的大小不符合规定。按照规则输入参数Window，详情请键入 help Median_Filter')
end
elseif strcmp(Window_type,'square')
if ~((size(Window,1) == 1 && size(Window,2) == 1) || (size(Window,1) == 3 && size(Window,2) == 1))
error('Window矩阵的大小不符合规定。按照规则输入参数Window，详情请键入 help Median_Filter')
end
else
if ~((size(Window,1) == 1 && size(Window,2) == 2) || (size(Window,1) == 3 && size(Window,2) == 2))
error('Window矩阵的大小不符合规定。按照规则输入参数Window，详情请键入 help Median_Filter')
end
end
% 规定窗口大小必须为奇数
if mean(mean(round(mod(Window,2)))) ~= 1
error('规定窗口的长宽必须为奇数')
end
% 请按规定输入参数Edge_Check。详情请键入 help Median_Filter
if ~strcmp(Edge_Check,'none')
error('请按规定输入参数Edge_Check。详情请键入 help Median_Filter')
end
if strcmp(Edge_Check,'none')
warning('没有对于函数或图像边缘的检测和处理！')
end
% 请按规定输入参数Drawing。详情请键入 help Median_Filter
if ~(strcmp(Drawing,'on') || strcmp(Drawing,'off'))
error('请按规定输入参数Drawing。详情请键入 help Median_Filter')
end
%% 函数核心部分
if Method == 1
if isempty(Input1)
Input1 = 1:length(Input2);
end
New = zeros(size(Input1,1),size(Input1,2));
arm = floor(Window/2);
area = arm+1:length(Input1)-arm;
if strcmp(Edge_Check,'none') % 无边缘处理
for i = 1:arm
New(1,i) = Input2(1,i);
end
for i = length(Input1)-arm+1:length(Input1)
New(1,i) = Input2(1,i);
end
end
%%%%%%%%%%%%%%%%%%%
% 此处后续可添加边缘处理
%%%%%%%%%%%%%%%%%%%
for i = area % 核心
M = Input2(1,i-arm:i+arm);
New(1,i) = median(M);
end
Output1 = Input1; % 输出
Output2 = New;
if strcmp(Drawing,'on') % 检测是否绘图
figure
subplot(121)
plot(Input1,Input2,'b-','LineWidth',1)
grid on
title('原输入')
subplot(122)
plot(Output1,Output2,'b-','LineWidth',1)
grid on
title('中值滤波后')
end
end
if Method == 2
New = zeros(size(Input1,1),size(Input1,2));
if strcmp(Window_type,'square') || strcmp(Window_type,'rectangle') || strcmp(Window_type,'cross')
if size(Window,2) == 1
r_arm = floor(Window/2);
c_arm = r_arm;
elseif size(Window,2) == 2
r_arm = floor(Window(1,2)/2); % 纵臂长
c_arm = floor(Window(1,1)/2); % 横臂长
end
end
r_area = r_arm+1:r_length-r_arm;
c_area = c_arm+1:c_length-c_arm;
if strcmp(Edge_Check,'none') % 无边缘处理
for i = 1:r_arm
New(i,:) = Input1(i,:);
end
for i = r_length-r_arm+1:r_length
New(i,:) = Input1(i,:);
end
for j = 1:c_arm
New(:,j) = Input1(:,j);
end
for j = c_length-c_arm+1:c_length
New(:,j) = Input1(:,j);
end
end
%%%%%%%%%%%%%%%%%%%
% 此处后续可添加边缘处理
%%%%%%%%%%%%%%%%%%%
if strcmp(Window_type,'square') || strcmp(Window_type,'rectangle') % 核心
for i = r_area
for j = c_area
M = Input1(i-r_arm:i+r_arm,j-c_arm:j+c_arm);
New(i,j) = median(median(M));
end
end
elseif strcmp(Window_type,'cross')
for i = r_area
for j = c_area
M = [Input1(i-r_arm:i+r_arm,j)',Input1(i,j-c_arm:j+c_arm)];
New(i,j) = median(median(M));
end
end
end
Output1 = New; % 输出
Output1 = uint8(Output1);
Output2 = [];
if strcmp(Drawing,'on') % 检测是否绘图
subplot(121)
imshow(Input1)
title('原图片')
subplot(122)
imshow(Output1)
title('中值滤波后图片')
end
end
if Method == 3
New = zeros(r_length,c_length,3);
if strcmp(Window_type,'square') || strcmp(Window_type,'rectangle') || strcmp(Window_type,'cross')
if size(Window,2) == 1
if size(Window,1) == 1
R_r_arm = floor(Window/2); % 红色纵臂长
R_c_arm = R_r_arm; % 红色横臂长
G_r_arm = floor(Window/2); % 绿色纵臂长
G_c_arm = G_r_arm; % 绿色横臂长
B_r_arm = floor(Window/2); % 蓝色纵臂长
B_c_arm = B_r_arm; % 蓝色横臂长
elseif size(Window,1) == 3
R_r_arm = floor(Window(1,1)/2); % 红色纵臂长
R_c_arm = R_r_arm; % 红色横臂长
G_r_arm = floor(Window(2,1)/2); % 绿色纵臂长
G_c_arm = G_r_arm; % 绿色横臂长
B_r_arm = floor(Window(3,1)/2); % 蓝色纵臂长
B_c_arm = B_r_arm; % 蓝色横臂长
end
elseif size(Window,2) == 2
if size(Window,1) == 1
R_r_arm = floor(Window(1,2)/2); % 红色纵臂长
R_c_arm = floor(Window(1,1)/2); % 红色横臂长
G_r_arm = floor(Window(1,2)/2); % 绿色纵臂长
G_c_arm = floor(Window(1,1)/2); % 绿色横臂长
B_r_arm = floor(Window(1,2)/2); % 蓝色纵臂长
B_c_arm = floor(Window(1,1)/2); % 蓝色横臂长
elseif size(Window,1) == 3
R_r_arm = floor(Window(1,2)/2); % 红色纵臂长
R_c_arm = floor(Window(1,1)/2); % 红色横臂长
G_r_arm = floor(Window(2,2)/2); % 绿色纵臂长
G_c_arm = floor(Window(2,1)/2); % 绿色横臂长
B_r_arm = floor(Window(3,2)/2); % 蓝色纵臂长
B_c_arm = floor(Window(3,1)/2); % 蓝色横臂长
end
end
end
R_r_area = R_r_arm+1:r_length-R_r_arm;
R_c_area = R_c_arm+1:c_length-R_c_arm;
G_r_area = G_r_arm+1:r_length-G_r_arm;
G_c_area = G_c_arm+1:c_length-G_c_arm;
B_r_area = B_r_arm+1:r_length-B_r_arm;
B_c_area = B_c_arm+1:c_length-B_c_arm;
if strcmp(Edge_Check,'none') % 无边缘处理
for k = 1:3
if k == 1
r_arm = R_r_arm;
c_arm = R_c_arm;
elseif k == 2
r_arm = G_r_arm;
c_arm = G_c_arm;
else
r_arm = B_r_arm;
c_arm = B_c_arm;
end
for i = 1:r_arm
New(i,:,k) = Input1(i,:,k);
end
for i = r_length-r_arm+1:r_length
New(i,:,k) = Input1(i,:,k);
end
for j = 1:c_arm
New(:,j,k) = Input1(:,j,k);
end
for j = c_length-c_arm+1:c_length
New(:,j,k) = Input1(:,j,k);
end
end
end
%%%%%%%%%%%%%%%%%%%
% 此处后续可添加边缘处理
%%%%%%%%%%%%%%%%%%%
for k = 1:3 % 核心
if k == 1
r_arm = R_r_arm;
c_arm = R_c_arm;
r_area = R_r_area;
c_area = R_c_area;
elseif k == 2
r_arm = G_r_arm;
c_arm = G_c_arm;
r_area = G_r_area;
c_area = G_c_area;
else
r_arm = B_r_arm;
c_arm = B_c_arm;
r_area = B_r_area;
c_area = B_c_area;
end
if strcmp(Window_type,'square') || strcmp(Window_type,'rectangle') % 核心
for i = r_area
for j = c_area
M = Input1(i-r_arm:i+r_arm,j-c_arm:j+c_arm,k);
New(i,j,k) = median(median(M));
end
end
elseif strcmp(Window_type,'cross')
for i = r_area
for j = c_area
M = [Input1(i-r_arm:i+r_arm,j,k)',Input1(i,j-c_arm:j+c_arm,k)];
New(i,j,k) = median(median(M));
end
end
end
Output1 = New; % 输出
Output1 = uint8(Output1);
Output2 = [];
if strcmp(Drawing,'on') % 检测是否绘图
subplot(121)
imshow(Input1)
title('原图片')
subplot(122)
imshow(Output1)
title('中值滤波后图片')
end
end
end
end
S1.2.1.3 中值滤波函数的测试
测试该函数的滤波效果。
1. 对于噪声密度较低的序列
对序列进行第一次中值滤波，滤波效果较好：
clc,clear all,close all
x = 0:0.01:2*pi;
y = sin(x);
noise = randi([1,length(x)],1,100);
y(1,noise) = y(1,noise)+randi([-500,500],1,length(noise))/1000;
[ans1,ans2] = Median_Filter(1,x,y,5,'default','none','on'); % 5*5正方形窗口中值滤波
对序列进行第二次中值滤波，滤波结果基本恢复原信号：
[ans1,ans2] = Median_Filter(1,ans1,ans2,7,'default','none','on'); % 7*7正方形窗口中值滤波
2. 对于噪声密度较高的序列
第一次滤波效果不显著：
clc,clear all,close all
x = 0:0.01:2*pi;
y = sin(x).*cos(2*x)+0.2*rand(1,length(x));
[ans1,ans2] = Median_Filter(1,x,y,3,'default','none','on');
第二次滤波后大体能恢复原信号形状，但是效果不佳：
[ans1,ans2] = Median_Filter(1,ans1,ans2,11,'default','none','on');
3. 对于有椒盐噪声的灰度图像
第一次滤波效果相当不错：
clc,clear all,close all
t = imread('D:\Huawei Share\OneHop\IMG_20210912_122924.jpg');
t = rgb2gray(t);
t = imnoise(t,'salt & pepper',0.1); % 加入密度为0.1的椒盐噪声
[ans1,ans2] = Median_Filter(2,t,[],3,'square','none','on'); % 3*3正方形窗口中值滤波
第二次滤波基本除去了大部分椒盐像素点：
[ans3,ans4] = Median_Filter(2,ans1,[],3,'square','none','on'); % 第二次中值滤波
4. 对于有椒盐噪声的RGB图像
图像滤波效果显著：
clc,clear all,close all
t = imread('D:\Huawei Share\OneHop\IMG_20210912_122924.jpg');
t = imnoise(t,'salt & pepper',0.1); % 加入密度为0.1的椒盐噪声
[ans1,ans2] = Median_Filter(3,t,[],[3;3;3],'square','none','on'); % 3*3正方形窗口中值滤波
代码由于没有经过优化，对于图像处理运行的时间会偏长。
S1.2.2 算术平均值滤波
算术平均值滤波只是在函数的核心段与中值滤波有区别。
function [Output1, Output2] = Mean_Filter(Method, Input1, Input2, Window, Window_type, Edge_Check, Drawing)
%
% Return the mean filtering result for binary image, RGB image
% or function image in rectangular coordinate system.
%
% Mean_Filter函数用于返回对于二值图像、RGB图像或直角坐标系下函数图像的算术平均值滤波结果
%
% 输入参数Method：取值可为1、2、3，其它值的输入会报错
% > Method = 1 时对直角坐标系下的函数图像进行算术平均值滤波
% > Method = 2 时对二值图像进行算术平均值滤波
% > Method = 3 时对RGB图像进行算术平均值滤波
% 输入参数Input：
% > 若Method = 1，则Input1、Input2分别输入函数对应的横纵坐标取值（均为行向量），
% 或者Input1输入空矩阵，Input2输入纵坐标取值，该函数将会从1开始自动填充横坐标
% > 若Method = 2，则Input1输入二值图像的二维矩阵，Input2输入空矩阵
% > 若Method = 3，则Input1输入RGB图像的三维矩阵
% 输入参数Window：算术平均值滤波的窗口大小，（其矩阵中的元素应均为正整数，若输入为双精度数则会向下取整）
% 输入参数Window_type：窗口类型
% > 若Window_type = 'default'（即窗口为一个行向量）
% ·Method = 1时，Window取一个小于Input2向量长度的整数，表示窗口长度（应为大于1的奇数）
% ·Method取其它值时会报错
% > 若Window_type = 'square'（即窗口为正方形）
% ·Method = 1时会报错
% ·若Method = 2或3，Window取一个小于Input1矩阵行宽的整数，表示窗口边长（应为大于1的奇数）
% > 若Window_type = 'rectangle'（即窗口为矩形）
% ·Method = 1时会报错
% ·若Method = 2或3，Window取一个小于Input1矩阵行宽的1*2向量，分别表示窗口长、宽（均应为大于1的奇数）
% > 若Window_type = 'cross'（即窗口为十字形）
% ·Method = 1时会报错
% ·若Method = 2或3，Window取一个小于Input1矩阵行宽的1*2向量，分别表示十字窗口长、宽（均应为大于1的奇数）
% > 若Method = 3，如果需要对R、G、B进行不同窗口大小的算术平均值滤波，请输入3*2或3*1的矩阵（视Window_type而定），
% 每一行分别表示对R、G、B进行算术平均值滤波处理的窗口大小
% 输入参数Edge_Check：边缘检测方式
% > 'none'（默认）：不进行边缘检测
% 输入参数Drawing：是否在函数内进行绘图或展示图像
% > 'off'（默认）：不展示图像
% > 'on'：展示图像
%
% 输出参数Output1：Method = 1时，输出处理结果的横坐标向量；Method = 2或3时输出图像对应的矩阵
% 输出参数Output2：Method = 1时，输出处理结果的纵坐标向量；Method = 2或3时Output2为空
%
% 常用调用方式：
% [Output1, Output2] = Mean_Filter(Method, Input1, Input2, Window, Window_type)
% [Output1, Output2] = Mean_Filter(Method, Input1, Input2, Window, Window_type, Edge_Check, Drawing)
%
% 编写于2021年10月6日20:32 程序员：K2SO4钾（华中师范大学）
%% 初判断错误与警告
% 设定初始参数
if isempty(Input1)
r_length = size(Input2,1);
c_length = size(Input2,2);
else
r_length = size(Input1,1);
c_length = size(Input1,2);
end
if nargin == 4
Window_type = 'default';
Edge_Check = 'none';
Drawing = 'off';
elseif nargin == 5
Edge_Check = 'none';
Drawing = 'off';
elseif nargin == 6
if strcmp(Edge_Check,'on') || strcmp(Edge_Check,'off')
Drawing = Edge_Check;
Edge_Check = 'none';
else
Drawing = 'off';
end
end
% 输入输出参数规范
if nargin >= 8
error('输入参数数量太多')
elseif nargin <= 3
error('输入参数数量太少')
end
if nargout >= 3
error('输出参数数量太多')
end
% Method参数规范
Method = floor(Method);
if Method >= 4 || Method < 1
error('Method的值规定为1、2、3，其它向下取整后在规定值之外的输入值将会报错')
end
% Method = 1时规定Input1、Input2应为等长的行向量
if Method == 1
if ~isempty(Input1)
if size(Input1,1) ~= 1 || size(Input2,1) ~= 1
error('Input1不为空时，Input1、Input2应为等长的行向量')
end
if length(Input1) ~= length(Input2)
error('Input1不为空时，Input1、Input2应为等长的行向量')
end
else
if size(Input2,1) ~= 1
error('Input2应为行向量')
end
end
end
% Method = 2、3时规定Input1的长宽必须在3*3以上
if Method == 2 || Method == 3
if size(Input1,1) <= 3 || size(Input1,2) <= 3
error('Method = 2、3时规定Input1的长宽必须在3*3以上')
end
if ~isempty(Input2)
warning('将忽略Input2的输入')
end
end
% 窗口大小应设定为正整数、大于1的奇数，且窗口长度或长宽应小于等于函数因变量长度或图像长宽
Window = floor(Window);
if sum(sum(Window < 3)) ~= 0
error('窗口大小应设定为正整数、大于1的奇数')
end
if Method == 1
if Window > c_length
error(error('窗口长度或长宽应小于等于函数因变量长度或图像长宽'))
end
elseif strcmp(Window_type,'square')
if max(Window(:,1)) > r_length || max(Window(:,1)) > c_length
error(error('窗口长度或长宽应小于等于函数因变量长度或图像长宽'))
end
else
if max(Window(:,1)) > r_length || max(Window(:,2)) > c_length
error('窗口长度或长宽应小于等于函数因变量长度或图像长宽')
end
end
% Method与Window_type搭配不符合规定。详情请键入 help Mean_Filter
if Method == 1
if ~(strcmp(Window_type,'default'))
error('Method与Window_type搭配不符合规定。详情请键入 help Mean_Filter')
end
else
if ~((strcmp(Window_type,'square')) || (strcmp(Window_type,'rectangle')) || (strcmp(Window_type,'cross')))
error('Method与Window_type搭配不符合规定。详情请键入 help Mean_Filter')
end
end
% Window矩阵的大小不符合规定。按照规则输入参数Window，详情请键入 help Mean_Filter
if Method == 1
if size(Window,1)*size(Window,2) ~= 1
error('Window矩阵的大小不符合规定。按照规则输入参数Window，详情请键入 help Mean_Filter')
end
elseif strcmp(Window_type,'square')
if ~((size(Window,1) == 1 && size(Window,2) == 1) || (size(Window,1) == 3 && size(Window,2) == 1))
error('Window矩阵的大小不符合规定。按照规则输入参数Window，详情请键入 help Mean_Filter')
end
else
if ~((size(Window,1) == 1 && size(Window,2) == 2) || (size(Window,1) == 3 && size(Window,2) == 2))
error('Window矩阵的大小不符合规定。按照规则输入参数Window，详情请键入 help Mean_Filter')
end
end
% 规定窗口大小必须为奇数
if mean(mean(round(mod(Window,2)))) ~= 1
error('规定窗口的长宽必须为奇数')
end
% 请按规定输入参数Edge_Check。详情请键入 help Mean_Filter
if ~strcmp(Edge_Check,'none')
error('请按规定输入参数Edge_Check。详情请键入 help Mean_Filter')
end
if strcmp(Edge_Check,'none')
warning('没有对于函数或图像边缘的检测和处理！')
end
% 请按规定输入参数Drawing。详情请键入 help Mean_Filter
if ~(strcmp(Drawing,'on') || strcmp(Drawing,'off'))
error('请按规定输入参数Drawing。详情请键入 help Mean_Filter')
end
%% 函数核心部分
if Method == 1
if isempty(Input1)
Input1 = 1:length(Input2);
end
New = zeros(size(Input1,1),size(Input1,2));
arm = floor(Window/2);
area = arm+1:length(Input1)-arm;
if strcmp(Edge_Check,'none') % 无边缘处理
for i = 1:arm
New(1,i) = Input2(1,i);
end
for i = length(Input1)-arm+1:length(Input1)
New(1,i) = Input2(1,i);
end
end
%%%%%%%%%%%%%%%%%%%
% 此处后续可添加边缘处理
%%%%%%%%%%%%%%%%%%%
for i = area % 核心
M = Input2(1,i-arm:i+arm);
New(1,i) = mean(M);
end
Output1 = Input1; % 输出
Output2 = New;
if strcmp(Drawing,'on') % 检测是否绘图
figure
subplot(121)
plot(Input1,Input2,'b-','LineWidth',1)
grid on
title('原输入')
subplot(122)
plot(Output1,Output2,'b-','LineWidth',1)
grid on
title('算术平均值滤波后')
end
end
if Method == 2
New = zeros(size(Input1,1),size(Input1,2));
if strcmp(Window_type,'square') || strcmp(Window_type,'rectangle') || strcmp(Window_type,'cross')
if size(Window,2) == 1
r_arm = floor(Window/2);
c_arm = r_arm;
elseif size(Window,2) == 2
r_arm = floor(Window(1,2)/2); % 纵臂长
c_arm = floor(Window(1,1)/2); % 横臂长
end
end
r_area = r_arm+1:r_length-r_arm;
c_area = c_arm+1:c_length-c_arm;
if strcmp(Edge_Check,'none') % 无边缘处理
for i = 1:r_arm
New(i,:) = Input1(i,:);
end
for i = r_length-r_arm+1:r_length
New(i,:) = Input1(i,:);
end
for j = 1:c_arm
New(:,j) = Input1(:,j);
end
for j = c_length-c_arm+1:c_length
New(:,j) = Input1(:,j);
end
end
%%%%%%%%%%%%%%%%%%%
% 此处后续可添加边缘处理
%%%%%%%%%%%%%%%%%%%
if strcmp(Window_type,'square') || strcmp(Window_type,'rectangle') % 核心
for i = r_area
for j = c_area
M = Input1(i-r_arm:i+r_arm,j-c_arm:j+c_arm);
New(i,j) = mean(mean(M));
end
end
elseif strcmp(Window_type,'cross')
for i = r_area
for j = c_area
M = [Input1(i-r_arm:i+r_arm,j)',Input1(i,j-c_arm:j+c_arm)];
New(i,j) = mean(mean(M));
end
end
end
Output1 = New; % 输出
Output1 = uint8(Output1);
Output2 = [];
if strcmp(Drawing,'on') % 检测是否绘图
subplot(121)
imshow(Input1)
title('原图片')
subplot(122)
imshow(Output1)
title('算术平均值滤波后图片')
end
end
if Method == 3
New = zeros(r_length,c_length,3);
if strcmp(Window_type,'square') || strcmp(Window_type,'rectangle') || strcmp(Window_type,'cross')
if size(Window,2) == 1
if size(Window,1) == 1
R_r_arm = floor(Window/2); % 红色纵臂长
R_c_arm = R_r_arm; % 红色横臂长
G_r_arm = floor(Window/2); % 绿色纵臂长
G_c_arm = G_r_arm; % 绿色横臂长
B_r_arm = floor(Window/2); % 蓝色纵臂长
B_c_arm = B_r_arm; % 蓝色横臂长
elseif size(Window,1) == 3
R_r_arm = floor(Window(1,1)/2); % 红色纵臂长
R_c_arm = R_r_arm; % 红色横臂长
G_r_arm = floor(Window(2,1)/2); % 绿色纵臂长
G_c_arm = G_r_arm; % 绿色横臂长
B_r_arm = floor(Window(3,1)/2); % 蓝色纵臂长
B_c_arm = B_r_arm; % 蓝色横臂长
end
elseif size(Window,2) == 2
if size(Window,1) == 1
R_r_arm = floor(Window(1,2)/2); % 红色纵臂长
R_c_arm = floor(Window(1,1)/2); % 红色横臂长
G_r_arm = floor(Window(1,2)/2); % 绿色纵臂长
G_c_arm = floor(Window(1,1)/2); % 绿色横臂长
B_r_arm = floor(Window(1,2)/2); % 蓝色纵臂长
B_c_arm = floor(Window(1,1)/2); % 蓝色横臂长
elseif size(Window,1) == 3
R_r_arm = floor(Window(1,2)/2); % 红色纵臂长
R_c_arm = floor(Window(1,1)/2); % 红色横臂长
G_r_arm = floor(Window(2,2)/2); % 绿色纵臂长
G_c_arm = floor(Window(2,1)/2); % 绿色横臂长
B_r_arm = floor(Window(3,2)/2); % 蓝色纵臂长
B_c_arm = floor(Window(3,1)/2); % 蓝色横臂长
end
end
end
R_r_area = R_r_arm+1:r_length-R_r_arm;
R_c_area = R_c_arm+1:c_length-R_c_arm;
G_r_area = G_r_arm+1:r_length-G_r_arm;
G_c_area = G_c_arm+1:c_length-G_c_arm;
B_r_area = B_r_arm+1:r_length-B_r_arm;
B_c_area = B_c_arm+1:c_length-B_c_arm;
if strcmp(Edge_Check,'none') % 无边缘处理
for k = 1:3
if k == 1
r_arm = R_r_arm;
c_arm = R_c_arm;
elseif k == 2
r_arm = G_r_arm;
c_arm = G_c_arm;
else
r_arm = B_r_arm;
c_arm = B_c_arm;
end
for i = 1:r_arm
New(i,:,k) = Input1(i,:,k);
end
for i = r_length-r_arm+1:r_length
New(i,:,k) = Input1(i,:,k);
end
for j = 1:c_arm
New(:,j,k) = Input1(:,j,k);
end
for j = c_length-c_arm+1:c_length
New(:,j,k) = Input1(:,j,k);
end
end
end
%%%%%%%%%%%%%%%%%%%
% 此处后续可添加边缘处理
%%%%%%%%%%%%%%%%%%%
for k = 1:3 % 核心
if k == 1
r_arm = R_r_arm;
c_arm = R_c_arm;
r_area = R_r_area;
c_area = R_c_area;
elseif k == 2
r_arm = G_r_arm;
c_arm = G_c_arm;
r_area = G_r_area;
c_area = G_c_area;
else
r_arm = B_r_arm;
c_arm = B_c_arm;
r_area = B_r_area;
c_area = B_c_area;
end
if strcmp(Window_type,'square') || strcmp(Window_type,'rectangle') % 核心
for i = r_area
for j = c_area
M = Input1(i-r_arm:i+r_arm,j-c_arm:j+c_arm,k);
New(i,j,k) = mean(mean(M));
end
end
elseif strcmp(Window_type,'cross')
for i = r_area
for j = c_area
M = [Input1(i-r_arm:i+r_arm,j,k)',Input1(i,j-c_arm:j+c_arm,k)];
New(i,j,k) = mean(mean(M));
end
end
end
Output1 = New; % 输出
Output1 = uint8(Output1);
Output2 = [];
if strcmp(Drawing,'on') % 检测是否绘图
subplot(121)
imshow(Input1)
title('原图片')
subplot(122)
imshow(Output1)
title('算术平均值滤波后图片')
end
end
end
end
clc,clear all,close all
t = imread('D:\Huawei Share\OneHop\IMG_20210912_122924.jpg');
t = imnoise(t,'salt & pepper',0.1); % 加入密度为0.1的椒盐噪声
[ans1,ans2] = Mean_Filter(3,t,[],[3,3;3,5;3,3],'rectangle','none','on'); % 3*3正方形窗口中值滤波
clc,clear all,close all
x = 0:0.01:2*pi;
y = sin(x).*cos(2*x)+0.2*rand(1,length(x));
[ans1,ans2] = Mean_Filter(1,x,y,9,'default','none','on');
[ans1,ans2] = Median_Filter(1,ans1,ans2,9,'default','none','on');
撰写：K2SO4钾（华中师范大学）
曹立学, 张鹏超. 计算机控制技术[M]. 西安：西安电子科技大学出版社. 2012. ↩︎
曹立学, 张鹏超. 计算机控制技术[M]. 西安：西安电子科技大学出版社. 2012. ↩︎
椒盐噪声 - 搜狗百科[EB/OL], doi: https://baike.sogou.com/v10913627.htm?fromTitle=%E6%A4%92%E7%9B%90%E5%99%AA%E5%A3%B0]
↩︎
扫一扫
9万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
