高斯高通滤波器对图像滤波

写在前面，这是一篇有关于"XX"大学光学"X"程的数值"XX"课作业的讲解汇总，江月年年望相似，作业也是如此。

本人中秋国庆回不了家，无事可干，故于2023年9月30日决定写下该回答，仅供参考，节约诸位时间。如果Anybody因为照抄导致你自己挂科，并且让老师抗拒这种形式的交流，那你是真该死

最后，愿本回答会成为帮助各位师弟、师妹的渡河之桥，更快抵达你的诗和远方

第一部分：实现的功能

A.要求功能一

1. 读入一幅灰度级不协调的图像，如pout.tif
   a.分析其直方图(文字说明)；
   b.设计灰度变换公式；
   c.调用imadjust函数，调整变换表达式参数到图像的灰度级分布均衡。
2. 对自选一幅图像调价噪声，使用均值滤波和中值滤波去噪，并显示和保存

B.要求功能二

1. 用傅立叶变换求由两个时间周期分别为0.2 和5的正弦函数;函数，y1，y2之和组成的函数的频谐。y1和y2的振幅分别为5和1。
   y1=5 * sin(2 * pi * 5 * x)-----y2 = sin(2 * pi * 0.2 * x)------y= y1+y2
2. a.在rice.tif中加噪声，并用2阶Butterworth低通滤波器去除噪声.
   b.用高斯高通滤波器对rice.tif原图图像滤波

第二部分：基本说明

A.用到的函数

imhist()

基本用法：用于显示直方图，要创建显示窗口

figure('Name','直方图');     %创建显示直方图的窗口
imhist(pic);                 %显示直方图

fspecial()

基本用法：fspecial(‘要的滤波器类型’,参数)

fil = fspecial('average',[3,3]);          %创建3 * 3的均值滤波矩阵。

滤波器字符	滤波器类型
‘average’	平均值滤波器
‘disk’	圆形平均值滤波器 (pillbox)
‘gaussian’	高斯低通滤波器。不推荐。请改用 imgaussfilt 或 imgaussfilt3。
‘laplacian’	逼近二维拉普拉斯算子
‘log’	高斯拉普拉斯滤波器
‘motion’	逼近相机的线性运动
‘prewitt’	Prewitt 水平边缘强调滤波器
‘sobel’	Sobel 水平边缘强调滤波器

imfilter()

基本用法：imfilter函数叫做实现线性空间滤波函数，功能是对任意类型数组或多维图像进行滤波，函数形式是 pic = imfilter(原图片矩阵,滤波矩阵)。

fil = fspecial('average',[3,3]);    %创建3 * 3的滤波矩阵。
pic_ave = imfilter(pic_noise,fil,'symmetric');   %均值滤波，边界之外的输入数组值是通过沿数组边界对数组进行镜面反射得到

imadjust()

基本用法：imadjust（f,[low_in high_in],[low_out high_out],gamma）

pic_ad = imadjust(pic,[0.29,0.63],[0,255]);     
%调整矩阵灰度范围，获得输入图片0.29*255到0.63*255的灰度值，线性扩展变为灰度值0到255的输出图

fft()、fft2()

基本用法：傅里叶变换函数，就是简单的直接用

y = [1,2,3,4,5,6,7,8,9,10]; %原函数,或者说一阶的行列式
Fy = fft(y);                %对原函数进行1阶快速傅里叶变换

pic = [1,2,3;3,2,1;6,6,6];  %2阶矩阵
FFT_pic = fft2(pic);        %对原图矩阵进行2阶离散傅里叶变换

结果展示：

%  y = [1,2,3,4,5,6,7,8,9,10];
Fy =

  21.0000 + 0.0000i  -3.0000 + 5.1962i  -3.0000 + 1.7321i  -3.0000 + 0.0000i  -3.0000 - 1.7321i  -3.0000 - 5.1962i

%   pic = [1,2,3;3,2,1;6,6,6];
FFT_pic =

  30.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
  -6.0000 +10.3923i  -3.0000 + 0.0000i  -1.5000 - 2.5981i
  -6.0000 -10.3923i  -1.5000 + 2.5981i  -3.0000 - 0.0000i

ifft()、ifft2()

基本用法：傅里叶逆变换函数，直接用就好了，同上

fftshift()、iffshift()

基本用法：它本质就是一个矩阵交换的方法，目的是将零频量移到中间。举个例子：

%ffftshift和ifftshift的例子
% 元素个数为偶数
x=[1,2,3,4,5,6];
x_fftshift=fftshift(x)
x_ifftshift=ifftshift(x)

% 元素个数为奇数
y=[1,2,3,4,5];
y_fftshift=fftshift(y)
y_ifftshift=ifftshift(y)

运行结果

%运行结果为
%########           x=[1,2,3,4,5,6];   元素个数为偶数
x_fftshift =

     4     5     6     1     2     3


x_ifftshift =

     4     5     6     1     2     3

%########           y=[1,2,3,4,5];    元素个数为奇数
y_fftshift =

     4     5     1     2     3


y_ifftshift =

     3     4     5     1     2

B.难点思路分析

首先，我觉得这东西没啥难点，和同学交流，发现他们对时间周期，采样间隔，采样频率之类的概念有巨多的疑问，所以我就先说说这个问题

1、概念理解

首先，你要知道对于计算机，没有常规上的变量而言，有的只是一个个有对应关系的值。对于一个sin(2 * pi * x)函数的例子，获得的永远是一个个有一定对应关系的值，每一个 t 的值确定下来，y 就确定了

x = 0:0.1:1
fx = sin(2*pi*x) 

figure;
plot(x,y,'r*');

%  运行结果展示：
x =

         0    0.1000    0.2000    0.3000    0.4000    0.5000    0.6000    0.7000    0.8000    0.9000    1.0000


fx =

         0    0.5878    0.9511    0.9511    0.5878    0.0000   -0.5878   -0.9511   -0.9511   -0.5878   -0.0000

也就是说，傅里叶变换中的 f(x) 其实是一个个数值，x 是从 0 到 N-1的一个个整数。带入一个 u 的值，将每一个 x 和 f(x) 累加起来，就会得到一个 F(u) 的值。
计算机中的傅里叶变换，要注意的是，f(x) 中的 x 对matlab来说，其实是数组的第几个元素。 也就是说 x 一定是从 1 到 N 。与你取的 x 值毫无关系

这是计算机进行傅里叶变换的基本逻辑，所以怎么取好 x 的值很重要，主要表现在对 fx 数组值的影响上。

A.采样间隔(核心)

对于上面sin(2 * pi * t)函数的例子，我们获得 ft 数组，对其傅里叶变换最重要的是， ft 数组中的值有没有涵盖sin函数中有特征的值，ft 数值的值。下图展示一个典型的采样间隔错误

对于sin(2 * pi * t)我们 t 的采样间隔至少要大于0.5，这样才能大概取得函数中有鲜明特征的值

采样频率 = 1 / 采样间隔；

采样数 = 周期长度 / 采样间隔 + 1；（加1 是加上开始的原点，即 100个间隔，要采101个点）

2、滤波器去噪声的实现

1、傅里叶变换，获得频谱
2、用滤波器作用频谱，即滤波器乘以频谱，获得滤波后频谱
3、对滤波后频谱进行傅里叶逆变换

A.循环

我交的作业就是用这个，因为这个会提高代码的可读性，你去看下我的作业第二题，你就知道了，就不重复思路代码了。

B.矩阵基本积(哈达玛积)

% 在rice.tif中加噪声，并用2阶Butterworth低通滤波器去除噪声. 
pic_path = 'rice.png' ;    %输入图片路径
pic = imread(pic_path);    %打开图片
pic=im2double(pic);        %图片类转为double类
figure();                  %创建显示窗口
subplot(1,2,1);            %组合绘图1
title('原图');             %原图标题
imshow(pic);               %显示图片
%#########       在rice.tif中加噪声      ###############
pic_noise = imnoise(pic,'salt & pepper');    %添加椒盐噪声。高斯噪声'gaussian'；泊松噪声'poisson'；高斯白噪声'localvar'；乘性噪声'speckle'。
subplot(1,2,2);             %组合绘图2
title('噪声图');            %噪声图标题
imshow(pic_noise);          %显示图片
%#########          傅里叶变换            ############
FFT_pic = fft2(pic);                  %原图矩阵进行2阶傅里叶变换
FFT_pic_shift = fftshift(fft2(pic));  %变换后矩阵进行频移
[u_max,v_max] = size(FFT_pic_shift);  %获取变换后矩阵的行列数

%#########   定义频谱中的参量u、v和中心频率   ##########
n = 2;                     %Butterworth的阶数
D0 = 20;                   %Butterworth的截止频率
u0 = u_max / 2;            %获取中心频率
v0 = v_max / 2;            %获取中心频率
u = 1 : u_max;             %定义 FFT_pic(u，v) 中的频率 u 的行列式
v = 1 : v_max;             %定义 FFT_pic(u，v) 中的频率 v 的行列式

%########       矩阵基本积思路         ##################
D=sqrt( (u - u0) .^ 2 + (v - v0) .^ 2);      %定义一个频域滤波的范围
H = 1 ./ (1 + (D ./ D0) .^ (2*n));           %定义一个对 频域滤波范围内外 的作用方法--即--滤波器
FFT_pic_fre = fftshift(FFT_pic_shift .* H);  %用滤波器作用 频移后的频谱 ，并再次进行频移

pic_butter = abs(ifft2(FFT_pic_fre));        %对矩阵进行逆变换后求绝对值

figure('Name','变换后butter降噪图');          %创建显示变换后降噪图的窗口
imshow(pic_butter);                          %显示图片

第三部分：我的作业

第一题

%   文件名：     DSH_2_1.m
%   版本：       1.0
%   创建时间：   2023年09月19日 星期二 18：37 
%   功能：       1. 读入一幅灰度级不协调的图像，如pout.tif
%                 a.分析其直方图(文字说明)；b.设计灰度变换公式；c.调用imadjust函数，调整变换表达式参数到图像的灰度级分布均衡。
%               2.对自选一幅图像调价噪声，使用均值滤波和中值滤波去噪，并显示和保存


%功能：     1. 读入一幅灰度级不协调的图像，如pout.tif
pic_path = 'pout.tif' ;    %输入图片路径

pic = imread(pic_path);    %打开图片

figure('Name','原图');     %创建显示原图的窗口
imshow(pic);               %显示图片

figure('Name','直方图');     %创建显示直方图的窗口
imhist(pic);                %显示直方图

%分析：直方图横坐标为灰度，纵坐标为灰度点个数，灰度大多数集中在75-175，所以我们要让其灰度范围变广，不那么集中。此外图片的灰度变化不明显，人眼不能很好的区分，要让各灰度间隔变大。


%设计公式：公式为线性变化，灰度范围由[75,160]线性变化到[0,255]，变换公式应该为   y = (255 / (160 - 75)) * (x - 75),化简为y = 3 * x -255 。

pic_ad = imadjust(pic,[0.29,0.63],[]);     %调整矩阵灰度范围
figure('Name','调整图');     %创建显示直方图的窗口
imshow(pic_ad);                 %显示直方图

figure('Name','调整直后方图');     %创建显示调整后直方图的窗口
imhist(pic_ad);                %显示直方图

%功能：    2.对自选一幅图像调价噪声，使用均值滤波和中值滤波去噪，并显示和保存
pic_path = 'pout.tif' ;    %输入图片路径
win_path = 'C:\Users\DSH\Desktop\';    %输入保存路径--默认是DSH的桌面
pic = imread(pic_path);    %打开图片

pic_noise = imnoise(pic,'salt & pepper');  %添加椒盐噪声。高斯噪声'gaussian'；泊松噪声'poisson'；高斯白噪声'localvar'；乘性噪声'speckle'。

figure('Name','噪声图');     %创建显示噪声图的窗口
imshow(pic_noise);               %显示图片

fil = fspecial('average',[3,3]);    %创建3 * 3的滤波矩阵。
pic_ave = imfilter(pic_noise,fil,'symmetric');   %均值滤波，边界之外的输入数组值是通过沿数组边界对数组进行镜面反射得到

figure('Name','均值降噪图');     %创建显示均值降噪图的窗口
imshow(pic_ave);               %显示图片
imwrite(pic_ave,strcat(win_path,'pic_ave.tif'));   %保存均值降噪图到桌面

pic_med = medfilt2(pic_noise,[3,3]);     %中值滤波，3 * 3范围选一个

figure('Name','中值降噪图');     %创建显示中值降噪图的窗口
imshow(pic_med);               %显示图片
imwrite(pic_med,strcat(win_path,'pic_med.tif'));   %保存中值降噪图到桌面

fprintf('您好，世界\t来自DSH的问候\n');

第二题

%   文件名：     DSH_2_2.m
%   版本：       1.0
%   创建时间：   2023年09月20日 星期二 16：28
%   功能：      1. 用傅立叶变换求由两个时间周期分别为0.2 和5的正弦函数;函数，y1，y2之和组成的函数的频谐。y1和y2的振幅分别为5和1。
%                         y1=5 * sin(2 * pi * 5 * x)        y2 = sin(2 * pi * 0.2 * x)       y= y1+y2
%              2. a.在rice.tif中加噪声，并用2阶Butterworth低通滤波器去除噪声.
%                 b.用高斯高通滤波器对rice.tif原图图像滤波


% %1、用傅立叶变换求由两个时间周期分别为0.2 和5的正弦函数;函数，y1，y2之和组成的函数的频谐。y1和y2的振幅分别为5和1。 y1=5 * sin(2 * pi * 5 * x)        y2 = sin(2 * pi * 0.2 * x)       y= y1+y2
Fs = 20;        %两函数频率分别为 5 和 0.2 ，采样频率应该大于 2 * 5 = 10 Hz; 这里取 20Hz
dx = 1 / Fs;    %采样间隔为采样频率的倒数，即1/20；
L = 201;        %采样数乘以采样间隔要大于函数一个周期(分别为 0.2 和 5)，即采样数(L)>最大周期(5)/采样间隔(1/20);即 L > 100 ，这里取201
df = Fs / (L - 1);     %频率间隔，即画出的频谱图的间隔
x = (0 : (L-1)) * dx;  %x的取值从0*dx到(L-1)*dx，即 0,1/20,2/20,3/20,

y1 = 5 * sin(2 * pi * 5 * x);   %函数y1
y2 = sin(2 * pi * 0.2 * x);     %函数y2
y = y1 + y2;           %两函数之和

fy =abs(fftshift(fft(y)));    %fft傅里叶变换;之后fftshift将零频分量移动到数组中心，重新排列傅里叶变换;最后取其绝对值。
fx = - Fs / 2 : df :Fs/2;     %频谱从-Fs/2到Fs/2，间隔为df

figure('name','频谱图');      %创建显示频谱图的窗口
plot(fx,fy);                  %显示频谱图

%2. a.在rice.tif中加噪声，并用2阶Butterworth低通滤波器去除噪声. b.用高斯高通滤波器对rice.tif原图图像滤波

pic_path = 'rice.png' ;    %输入图片路径
pic = imread(pic_path);    %打开图片
pic=im2double(pic);        %
figure();                  %创建显示窗口
subplot(1,2,1);            %组合绘图1
title('原图');             %原图标题
imshow(pic);               %显示图片
pic_noise = imnoise(pic,'salt & pepper');    %添加椒盐噪声。高斯噪声'gaussian'；泊松噪声'poisson'；高斯白噪声'localvar'；乘性噪声'speckle'。
subplot(1,2,2);             %组合绘图2
title('噪声图');            %噪声图标题
imshow(pic_noise);          %显示图片

FFT_pic = fft2(pic);        %原图矩阵进行2阶傅里叶变换
FFT_pic_shift = fftshift(fft2(pic));  %变换后矩阵进行频移

[u_max,v_max] = size(FFT_pic_shift);   %获取变换后矩阵的行列数

%#########a.在rice.tif中加噪声，并用2阶Butterworth低通滤波器去除噪声######################

n = 2;                     %Butterworth的阶数
D0 = 20;                   %Butterworth的截止频率
u0 = u_max / 2;            %获取中心频率
v0 = v_max / 2;            %获取中心频率
FFT_pic_butter_shift = ones(u_max,v_max);    %创建全一矩阵
for x=1:u_max
    for y=1:v_max
      
        D=sqrt((x-u0)^2+(y-v0)^2);  %计算点（x,y）到中心点的距离
      
        h=1/(1+(D/D0)^(2*n));       %计算巴特沃斯滤波器
    
        FFT_pic_butter_shift(x,y)=FFT_pic_shift(x,y)*h;    %用滤波器乘以主函数
    end
end
FFT_pic_butter = fftshift(FFT_pic_butter_shift); %对矩阵进行频移
pic_butter = abs(ifft2(FFT_pic_butter));    %对矩阵进行逆变换后求绝对值
figure('Name','变换后butter降噪图');        %创建显示变换后降噪图的窗口
imshow(pic_butter);                     %显示图片

%#########b.用高斯高通滤波器对rice.tif原图图像滤波#################################

u0 = u_max / 2;      %获取中心频率
v0 = v_max / 2;      %获取中心频率
D0 = 20;             %gauss的截止频率
FFT_pic_gauss_shift = ones(u_max,v_max);    %创建全一矩阵
for x=1:u_max
    for y=1:v_max
      
        D=sqrt((x-u0)^2+(y-v0)^2);  %计算点（x,y）到中心点的距离
      
        h = 1 - exp(-D^2 / (2 * D0 ^ 2));       %计算高斯滤波器
    
        FFT_pic_gauss_shift(x,y)=FFT_pic_shift(x,y)*h;    %用滤波器乘以主函数
    end
end
FFT_pic_gauss = fftshift(FFT_pic_gauss_shift); %对矩阵进行频移
pic_gauss = abs(ifft2(FFT_pic_gauss));    %对矩阵进行逆变换后求绝对值
figure('Name','变换后gauss降噪图');        %创建显示变换后降噪图的窗口
imshow(pic_gauss);                        %显示图片

fprintf('您好，世界\t来自DSH的问候\n');