数字图像处理实验（开源）

@HNUSTer

已于 2024-05-20 16:47:41 修改

阅读量981

点赞数 25

分类专栏：《数字图像处理（MATLAB版）》测绘工程专业课笔记/期末复习资料文章标签：数字图像处理实验图像读取图像傅里叶变换图像高通滤波图像低通滤波

于 2024-04-23 19:57:34 首次发布

本文链接：https://blog.csdn.net/zk180531/article/details/138135379

版权

测绘工程专业课笔记/期末复习资料同时被 2 个专栏收录

30 篇文章 66 订阅

订阅专栏

《数字图像处理（MATLAB版）》

5 篇文章 1 订阅

订阅专栏

1 图像读取

1.1 程序代码

clc
clear
close all
%%
img = imread('RedBull.jpg');
imshow(img);
title('运行结果');

1.2 运行结果

1.3 代码分析

clc: 这个命令用于清除MATLAB命令窗口的内容，以便于清除之前的输出，让新的输出更清晰。

clear: 这个命令用于清除工作区中的所有变量，以确保开始一个干净的工作环境。

close all: 这个命令关闭所有已经打开的图形窗口，以确保开始时没有任何窗口打开。

%%: 在MATLAB中，%% 符号用于标记代码节（section），它可以用来将代码分成不同的部分，提高可读性，并且可以方便地运行或调试部分代码。

img = imread('RedBull.jpg');: 这一行代码将名为 "RedBull.jpg" 的图像文件读取到名为 img 的变量中。imread 函数用于读取图像文件，它的参数是图像文件的路径，可以是相对路径或绝对路径。

imshow(img);: 这一行代码用于在图像窗口中显示 img 变量中存储的图像。imshow 函数用于显示图像，它的参数是要显示的图像数据。

title('运行结果');: 这一行代码用于为图像窗口添加标题，标题内容为 "运行结果"。

2 图像傅里叶变换

2.1 傅里叶级数可视化

2.1.1 程序代码

clc
clear
close all
%%
m = 15;   % Series
N = 512;  % Period
T0 = 4; 
T = 2;
% Time vectors
t1 = linspace(-T/2, T/2, N);
t2 = linspace(T/2, T0-T/2, N);
t  = [(t1-T0)'; (t2-T0)'; t1'; t2'; (t1+T0)'];
% Construct periodic signals
% s = zeros(5*N,1);
% s(1:N)=1; s(2*N+1:3*N)=1; s(4*N+1:5*N)=1;
s = [ones(N,1);zeros(N,1);ones(N,1);zeros(N,1);ones(N,1)]'; 
y = zeros(m, length(s));
ys = s-0.5;
hold on; grid on;
% Draw period square wave
plot3(t-(t+2), t, ys,'b','LineWidth',2.5); 
% Draw amplitude-frequency curve
A = 0.5;  % Coefficient
fsamp = 1024;
f = linspace(1, m+1, fsamp);
freq = 1:m;
F = A*sinc(A*f)*5;
plot3(f,T0+1+f*0,F,'r-','LineWidth',2.5)
fill3([f(1),f,f(end)],[T0+1,T0+1+f*0,T0+1],[min(F(:)),F,min(F(:))],...
    'r','FaceAlpha',0.2)
% Draw harmonic signals
x = A*zeros(size(t)); % Harmonic signals
X = A*ones(size(t));  % Synthesis signals
for k = 1:m
    x = 2*A*sinc(A*k)*cos(2*pi*t*k/T0);
    y(k,:) = x; 
    plot3(k+t*0,t,y(k,:),'LineWidth',1.5); % Harmonic
    X = X + x;
end
% Draw synthesis signals
plot3(k+1+t*0,t,X-0.5,'LineWidth',2.5,'Color','g');
xlabel('Frequency view','Rotation', 15);
ylabel('Time duration','Rotation', -10);
zlabel('Amplitude/Magnitude');
axis([-2,m+1,-(T0+1),T0+1,-1,2])
set(gca,'XTick',-T0-1:2:m);
set(gca,'YTick',-T0-1:1:T0+1);
set(gca,'ZTick',-1:0.5:T0);
title('傅里叶级数展开3D可视化','FontSize',15);
view(-49,23)
set(gcf,'Color','White');

2.1.2 运行结果

2.1.3 代码分析

初始化环境：
clc, clear, close all：清除命令窗口、清除工作区、关闭所有图形窗口，确保开始时环境清洁。
m = 15：系列数量，即要展示的傅里叶级数的级数。
N = 512：周期的采样点数。
T0 = 4：总周期。
T = 2：一个周期内的时间长度。

生成时间向量：
t1 和 t2 是两个分段时间向量，分别表示了一个周期内的两段时间。
t 是按照特定顺序连接了这些时间向量，以模拟周期性信号。

构造周期信号：
s 是一个周期信号，其中间段为 1，其余段为 0。
y 是一个矩阵，用于存储不同频率下的傅里叶级数展开结果。

绘制周期方波：
使用 plot3 函数绘制了周期方波图形，其中 t-(t+2) 是为了让图形在 x 轴上有所错位。

绘制幅频特性曲线：
计算并绘制了幅频特性曲线，即幅度随频率变化的曲线。
使用 fill3 函数填充了幅频特性曲线下的区域，使其更加直观。

绘制谐波信号：
使用循环计算并绘制了每个谐波信号。
sinc 函数生成了谐波信号的幅度。
使用 cos 函数生成了正弦波。
最后将每个谐波信号相加得到综合信号。

绘制合成信号：
绘制了所有谐波信号相加后的合成信号。

设置图形属性：
设置了坐标轴的标签和范围，以及标题。
设置了坐标轴的刻度。
调整了视角。

展示结果：
显示了绘制好的 3D 图形。

2.2 图像傅里叶变换

2.2.1 程序代码

clc
clear
close all
%%
img=imread('MTB01.png');
subplot(2,2,1);
imshow(img);
title('原始图像');
f=rgb2gray(img);    %对于RGB图像必须做的一步，也可以用im2double函数
F=fft2(f);          %傅里叶变换
F1=log(abs(F)+1);   %取模并进行缩放
subplot(2,2,2);
imshow(F1,[]);
title('傅里叶变换频谱图像');
Fs=fftshift(F);      %将频谱图中零频率成分移动至频谱图中心
S=log(abs(Fs)+1);    %取模并进行缩放
subplot(2,2,3);
imshow(S,[]);
title('频移后的频谱图像');
fr=real(ifft2(ifftshift(Fs)));  %频率域反变换到空间域，并取实部
ret=im2uint8(mat2gray(fr));    %更改图像类型
subplot(2,2,4);
imshow(ret),title('逆傅里叶变换图像');

2.2.2 运行结果

2.2.3 代码分析

clc, clear, close all: 这些是清除 MATLAB 命令窗口、清除工作区变量、关闭所有图形窗口的命令，确保开始时环境干净。

img=imread('MTB01.png');: 这一行代码读取名为 "MTB01.png" 的图像文件，并将其存储在名为 img 的变量中。

subplot(2,2,1); imshow(img); title('原始图像');: 这一部分创建了一个 2x2 的子图，并在第一个子图中显示了原始图像 img，并添加了标题 "原始图像"。

f=rgb2gray(img);: 如果输入图像是 RGB 彩色图像，这一行代码将其转换为灰度图像，存储在变量 f 中。

F=fft2(f);: 这一行代码对灰度图像进行二维傅里叶变换，结果存储在变量 F 中。

F1=log(abs(F)+1);: 这一行代码计算傅里叶变换结果的幅度谱，并对其取对数，以便更好地可视化。+1 是为了避免对数运算中出现零值。

subplot(2,2,2); imshow(F1,[]); title('傅里叶变换频谱图像');: 这一部分在第二个子图中显示了傅里叶变换后的频谱图像 F1，并添加了标题 "傅里叶变换频谱图像"。

Fs=fftshift(F);: 这一行代码将傅里叶变换结果中的零频率分量移动到频谱的中心。

S=log(abs(Fs)+1);: 这一行代码计算频移后的频谱图像，并对其取对数。

subplot(2,2,3); imshow(S,[]); title('频移后的频谱图像');: 这一部分在第三个子图中显示了频移后的频谱图像 S，并添加了标题 "频移后的频谱图像"。

fr=real(ifft2(ifftshift(Fs)));: 这一行代码对频移后的频谱进行逆傅里叶变换，然后取实部得到恢复的图像。

ret=im2uint8(mat2gray(fr));: 这一行代码将逆傅里叶变换得到的图像进行类型转换和灰度缩放，确保其值在 [0,255] 的范围内。

subplot(2,2,4); imshow(ret),title('逆傅里叶变换图像');: 这一部分在第四个子图中显示了逆傅里叶变换后的图像 ret，并添加了标题 "逆傅里叶变换图像"。

3 图像高通滤波

3.1 算法1——理想高通滤波

3.1.1 程序代码

clc
clear
close all
%%
%理想高通
I = imread('house.tif');
figure(1);
subplot(221),imshow(I);
title('原始图像');
I=im2double(I);
s=fftshift(fft2(I));%傅里叶变换，直流分量搬移到频谱中心
subplot(222), imshow(log(abs(s)+1),[]); 
title('傅里叶变换取对数所得频谱图像');
[a,b]=size(s);
h=zeros(a,b);%滤波器函数
res=zeros(a,b);%保存结果
a0=round(a/2);
b0=round(b/2);
d=40;
for i=1:a 
    for j=1:b 
        distance=sqrt((i-a0)^2+(j-b0)^2);
        if distance<d
            h(i,j)=0;
        else
            h(i,j)=1;
        end
    end
end
res=s.*h;
res=real(ifft2(ifftshift(res)));
subplot(223),imshow(res);
title('理想高通滤波所得图像'); 
subplot(224),imshow(h);
title('理想高通滤波器图像');

3.1.2 运行结果

3.1.3 代码分析

clc, clear, close all: 这些命令用于清除命令窗口、工作区和关闭所有图形窗口，以确保开始时处于一个干净的状态。

I = imread('house.tif');: 从名为 "house.tif" 的图像文件中读取图像数据，并将其存储在变量 I 中。

figure(1);: 创建一个新的图形窗口，图形窗口的编号为 1。

subplot(221),imshow(I);: 将当前图形窗口分割为 2 行 2 列，选择第一个子图，并在该子图中显示原始图像 I。

%理想高通: 这是一行注释，用于标识以下代码段的功能，表示接下来要实现理想高通滤波器。

I=im2double(I);: 将原始图像 I 转换为双精度浮点数格式。

s=fftshift(fft2(I));: 对双精度浮点数格式的图像 I 进行二维傅里叶变换，并通过 fftshift 函数将频谱的直流分量移到频谱中心，结果保存在变量 s 中。

subplot(222), imshow(log(abs(s)+1),[]);: 选择第二个子图，并在该子图中显示 s 的幅度谱，使用 log(abs(s)+1) 对幅度进行对数变换，以便更好地显示频谱图像。

[a,b]=size(s);: 获取频谱 s 的尺寸，并将其分别存储在变量 a 和 b 中。

h=zeros(a,b);: 创建一个大小与频谱 s 相同的零矩阵，用于存储理想高通滤波器的频率响应。

res=zeros(a,b);: 创建一个大小与频谱 s 相同的零矩阵，用于存储滤波后的频谱。

a0=round(a/2); b0=round(b/2); d=40;: 计算频谱中心坐标 (a0, b0)，以及理想高通滤波器的截止频率 d。

for i=1:a: 外层循环遍历频谱 s 的行数。

for j=1:b: 内层循环遍历频谱 s 的列数。

distance=sqrt((i-a0)^2+(j-b0)^2);: 计算当前频谱点到中心点的距离。

if distance<d: 如果距离小于截止频率 d，则将对应位置的滤波器系数设为 0，表示该频率成分被抑制。

else: 如果距离大于等于截止频率 d，则将对应位置的滤波器系数设为 1，表示该频率成分保留。

res=s.*h;: 将原始频谱 s 与滤波器 h 相乘，得到滤波后的频谱 res。

res=real(ifft2(ifftshift(res)));: 对滤波后的频谱 res 进行逆傅里叶变换，并通过 ifftshift 和 real 函数将结果转换为实数类型，得到滤波后的图像。

subplot(223),imshow(res);: 选择第三个子图，并在该子图中显示滤波后的图像。

title('理想高通滤波所得图像');: 设置第三个子图的标题。

subplot(224),imshow(h);: 选择第四个子图，并在该子图中显示理想高通滤波器的频率响应。

title('理想高通滤波器图像');: 设置第四个子图的标题。

3.2 算法2——高斯高通滤波

3.2.1 程序代码

clc
clear
close all
%%
%高斯高通
I=imread('a.tif');
subplot(221),imshow(I);
title('原始图像');
Y=fft2(im2double(I));%傅里叶变换
Y=fftshift(Y);%频谱搬移,直流分量搬移到频谱中心
subplot(222), imshow(log(abs(Y)+1),[]); 
title('傅里叶变换取对数所得频谱图像');
[M,N]=size(Y);%获得图像的高度和宽度
h=zeros(M,N);%滤波器函数
%图像中心点
M0=M/2;
N0=N/2;
%截至频率距离圆点的距离，delta表示高斯曲线的扩散程度
D0=40;
delta=D0;
for x=1:M
    for y=1:N
        %计算点（x,y）到中心点的距离
        d2=(x-M0)^2+(y-N0)^2;
        %计算高斯滤波器
        h(x,y)=1-exp(-d2/(2*delta^2));
    end
end
%滤波后结果
res=h.*Y;
res=real(ifft2(ifftshift(res)));
subplot(223),imshow(res);
title('高斯高通滤波所得图像'); 
subplot(224),imshow(h);
title('高斯高通滤波器图象');

3.2.2 运行结果

3.2.3 代码分析

clc、clear、close all：清除命令行窗口内容，清除工作区变量，关闭所有图形窗口。

I = imread('a.tif');：从名为 "a.tif" 的文件中读取图像数据，并将其存储在变量 I 中。这个图像文件应该在当前工作目录下。

subplot(221), imshow(I); title('原始图像');：将图像显示在 2x2 的子图中的第一个位置，并设置标题为 "原始图像"。

Y = fft2(im2double(I));：对图像 I 进行二维傅里叶变换，使用 im2double 函数将图像转换为双精度类型，结果存储在变量 Y 中。

Y = fftshift(Y);：将傅里叶变换后的频谱图进行频谱搬移，即将直流分量移到频谱中心。

subplot(222), imshow(log(abs(Y) + 1), []); title('傅里叶变换取对数所得频谱图像');：将经过对数变换的傅里叶频谱图像显示在子图中的第二个位置，并设置标题为 "傅里叶变换取对数所得频谱图像"。log(abs(Y) + 1) 是为了避免对数计算时出现对数值为零的情况，加上一个常数 1。

[M, N] = size(Y);：获取频谱图像的尺寸，其中 M 是行数，N 是列数。

h = zeros(M, N);：创建一个与频谱图像大小相同的全零矩阵，用于存储高斯滤波器函数。

M0 = M / 2; N0 = N / 2;：计算频谱图像中心点的坐标。

D0 = 40; delta = D0;：设置高斯滤波器的截止频率距离圆点的距离，delta 表示高斯曲线的扩散程度。

for x = 1:M 和 for y = 1:N：遍历频谱图像的所有像素点。

d2 = (x - M0)^2 + (y - N0)^2;：计算每个像素点到图像中心点的距离的平方。

h(x, y) = 1 - exp(-d2 / (2 * delta^2));：根据高斯函数的定义计算高斯滤波器的值，这里使用了二维高斯函数。

res = h .* Y;：将高斯滤波器与频谱图像相乘，得到滤波后的结果。

res = real(ifft2(ifftshift(res)));：对滤波后的结果进行逆傅里叶变换，并取实部，以获取空域图像。

subplot(223), imshow(res); title('高斯高通滤波所得图像');：将经过高斯高通滤波后的图像显示在子图中的第三个位置，并设置标题为 "高斯高通滤波所得图像"。

subplot(224), imshow(h); title('高斯高通滤波器图象');：将高斯高通滤波器的图像显示在子图中的第四个位置，并设置标题为 "高斯高通滤波器图象"。

3.3 算法3——巴特沃斯高通滤波

3.3.1 程序代码

clc
clear
close all
%%
%巴特沃斯高通
figure(3);
I=imread('fingerprint.tif');
subplot(221),imshow(I);
title('原始图像');
Y=fft2(im2double(I));%傅里叶变换
Y=fftshift(Y);%频谱搬移,直流分量搬移到频谱中心
subplot(222), imshow(log(abs(Y)+1),[]); 
title('傅里叶变换取对数所得频谱图像');
[M,N]=size(Y);%获得图像的高度和宽度
h=zeros(M,N);%滤波器函数
%图像中心点
M0=M/2;
N0=N/2;
d0=40;
%巴特沃斯滤波器的阶数
n_0=2;
for x=1:M
    for y=1:N
        distance=sqrt((x-M0)^2+(y-N0)^2);
        h(x,y)=1/(1+(d0/distance)^(2*n_0));
    end
end
%滤波后结果
res=h.*Y;
res=real(ifft2(ifftshift(res)));
subplot(223),imshow(res);
title('巴特沃斯高通滤波所得图像'); 
subplot(224),imshow(h);
title("巴特沃斯高通滤波器图象");

3.3.2 运行结果

3.3.3 代码分析

clc, clear, close all: 这三行分别用于清除命令窗口的内容、清除工作区的所有变量，关闭所有图形窗口。

%%: 在 MATLAB 中，%% 是用来标记代码段的分隔符，可以帮助将代码按照功能块分隔开来。

figure(3): 创建一个新的图形窗口，窗口编号为 3。

I=imread('fingerprint.tif');: 读取名为 fingerprint.tif 的指纹图像，并将其存储在变量 I 中。

subplot(221),imshow(I);: 将图像显示在窗口中，subplot(221) 表示将图像显示在 2x2 窗格中的第一个位置。

title('原始图像');: 设置窗口的标题为 "原始图像"。

Y=fft2(im2double(I));: 对原始图像进行二维傅里叶变换，并将结果存储在变量 Y 中。im2double 函数用于将图像转换为 double 类型，因为傅里叶变换要求输入是双精度浮点数。

Y=fftshift(Y);: 将频谱移动，使得频谱的中心对齐在图像的中心。

subplot(222), imshow(log(abs(Y)+1),[]);: 将经过傅里叶变换后的频谱图像显示在窗口中，log(abs(Y)+1) 用于增强频谱图像的对比度，[] 表示将图像的像素值映射到显示范围。

title('傅里叶变换取对数所得频谱图像');: 设置窗口的标题为 "傅里叶变换取对数所得频谱图像"。

[M,N]=size(Y);: 获取频谱图像的大小，M 表示行数，N 表示列数。

h=zeros(M,N);: 创建一个大小为 M x N 的全零矩阵，用于存储滤波器函数。

M0=M/2; N0=N/2;: 计算频谱图像的中心点坐标。

d0=40;: 设定巴特沃斯滤波器的截止频率。

n_0=2;: 设定巴特沃斯滤波器的阶数。

for x=1:M for y=1:N ... end end: 使用两层循环遍历频谱图像中的每一个像素点。

distance=sqrt((x-M0)^2+(y-N0)^2);: 计算当前像素点到频谱图像中心的距离。

h(x,y)=1/(1+(d0/distance)^(2*n_0));: 根据巴特沃斯高通滤波器的公式，计算滤波器函数的值。

res=h.*Y;: 将滤波器函数和频谱图像相乘，得到滤波后的结果。

res=real(ifft2(ifftshift(res)));: 对滤波后的结果进行逆傅里叶变换，然后取实部，得到最终的图像。

subplot(223),imshow(res);: 将滤波后的图像显示在窗口中。

title('巴特沃斯高通滤波所得图像');: 设置窗口的标题为 "巴特沃斯高通滤波所得图像"。

subplot(224),imshow(h);: 将巴特沃斯高通滤波器的图像显示在窗口中。

title("巴特沃斯高通滤波器图象");: 设置窗口的标题为 "巴特沃斯高通滤波器图象"。

4 图像低通滤波

4.1 算法1——理想低通滤波

4.1.1 程序代码

clc
clear
close all
%%
%理想低通
I = imread('house.tif');
figure(1);
subplot(221),imshow(I);
title('原始图像');
I=im2double(I);
s=fftshift(fft2(I));%傅里叶变换，直流分量搬移到频谱中心
subplot(222), imshow(log(abs(s)+1),[]); 
title('傅里叶变换取对数所得频谱图像');
[a,b]=size(s);
h=zeros(a,b);%滤波器函数
res=zeros(a,b);%保存结果
a0=round(a/2);
b0=round(b/2);
d=40;
for i=1:a 
    for j=1:b 
        distance=sqrt((i-a0)^2+(j-b0)^2);
        if distance<=d
            h(i,j)=1;
        else
            h(i,j)=0;
        end
    end
end
res=s.*h;
res=real(ifft2(ifftshift(res)));
subplot(223),imshow(res);
title('理想低通滤波所得图像'); 
subplot(224),imshow(h);
title("理想低通滤波器图象");

4.1.2 运行结果

4.1.3 代码分析

clc、clear、close all：这些命令用于清除 MATLAB 命令窗口的内容、清除工作区的所有变量和关闭所有图形窗口。

%：表示注释，后面的内容不会被 MATLAB 执行，只是用于对代码进行注解说明。

I = imread('house.tif');：从名为 'house.tif' 的图像文件中读取图像数据，并将其存储在变量 I 中。

figure(1);：创建一个新的图形窗口，并将其编号为 1。

subplot(221),imshow(I);：在当前图形窗口中创建一个 2x2 的子图，并将其设为第一个子图，然后显示原始图像 I。

title('原始图像');：给当前子图添加标题为 "原始图像"。

I=im2double(I);：将图像 I 转换为双精度浮点数表示，这通常是进行图像处理的一种常见做法。

s=fftshift(fft2(I));：对图像 I 进行二维傅里叶变换，并通过 fftshift 函数将结果进行频域中心化，将直流分量移动到频谱的中心。

subplot(222), imshow(log(abs(s)+1),[]);：在当前图形窗口中创建第二个子图，并显示傅里叶变换结果的幅度谱，取对数以增强显示效果。

title('傅里叶变换取对数所得频谱图像');：给当前子图添加标题。

[a,b]=size(s);：获取傅里叶变换结果 s 的大小，并将其分别存储在变量 a 和 b 中。

h=zeros(a,b);：创建一个与傅里叶变换结果相同大小的零矩阵，用于存储滤波器函数。

res=zeros(a,b);：创建一个与傅里叶变换结果相同大小的零矩阵，用于存储滤波后的结果。

a0=round(a/2); b0=round(b/2);：计算频谱中心的坐标。

d=40;：定义理想低通滤波器的截止频率。

for i=1:a for j=1:b：嵌套循环遍历傅里叶变换结果的每个像素。

distance=sqrt((i-a0)^2+(j-b0)^2);：计算当前像素到频谱中心的距离。

if distance<=d h(i,j)=1; else h(i,j)=0; end：根据距离与截止频率的比较，设定理想低通滤波器的滤波函数。

res=s.*h;：将原始频谱与滤波器函数进行点乘，实现频域滤波。

res=real(ifft2(ifftshift(res)));：对滤波后的结果进行逆傅里叶变换，并通过 real 函数取实部，得到滤波后的空域图像。

subplot(223),imshow(res);：在当前图形窗口中创建第三个子图，并显示滤波后的图像。

title('理想低通滤波所得图像');：给当前子图添加标题。

subplot(224),imshow(h);：在当前图形窗口中创建第四个子图，并显示理想低通滤波器的图像。

title("理想低通滤波器图象");：给当前子图添加标题。

4.2 算法2——高斯低通滤波

4.2.1 程序代码

clc
clear
close all
%%
%高斯低通
I=imread('a.tif');
subplot(221),imshow(I);
title('原始图像');
Y=fft2(im2double(I));%傅里叶变换
Y=fftshift(Y);%频谱搬移,直流分量搬移到频谱中心
subplot(222), imshow(log(abs(Y)+1),[]); 
title('傅里叶变换取对数所得频谱图像');
[M,N]=size(Y);%获得图像的高度和宽度
h=zeros(M,N);%滤波器函数
%图像中心点
M0=M/2;
N0=N/2;
%截至频率距离圆点的距离，delta表示高斯曲线的扩散程度
D0=40;
delta=D0;
for x=1:M
    for y=1:N
        %计算点（x,y）到中心点的距离
        d2=(x-M0)^2+(y-N0)^2;
        %计算高斯滤波器
        h(x,y)=exp(-d2/(2*delta^2));
    end
end
%滤波后结果
res=h.*Y;
res=real(ifft2(ifftshift(res)));
subplot(223),imshow(res);
title('高斯低通滤波所得图像'); 
subplot(224),imshow(h);
title("高斯低通滤波器图象");

4.2.2 运行结果

4.2.3 代码分析

clc, clear, close all: 清除命令窗口、清除工作区、关闭所有图像窗口，确保每次运行代码时都从一个干净的状态开始。

I=imread('a.tif');: 从名为"a.tif"的文件中读取图像，并将其存储在变量I中。

subplot(221),imshow(I); title('原始图像');: 在一个 2x2 的图像网格中创建第一个子图，并在该子图中显示原始图像I，同时设置标题为"原始图像"。

Y=fft2(im2double(I));: 对原始图像进行二维离散傅里叶变换（DFT），并将结果存储在变量Y中。im2double函数将图像的数据类型转换为双精度浮点数。

Y=fftshift(Y);: 对频谱进行移动，使得直流分量位于频谱的中心位置。

subplot(222), imshow(log(abs(Y)+1),[]); title('傅里叶变换取对数所得频谱图像');: 在图像网格中创建第二个子图，显示对数变换后的傅里叶频谱图像，log(abs(Y)+1)是为了避免对数运算时出现对数值为0的情况，[]表示使用全局图像灰度范围，title设置标题为"傅里叶变换取对数所得频谱图像"。

[M,N]=size(Y);: 获取傅里叶变换后频谱的大小，其中M表示行数，N表示列数。

h=zeros(M,N);: 创建一个大小为MxN的全零矩阵，用于存储高斯滤波器函数。

M0=M/2; N0=N/2;: 计算频谱中心点的坐标，M0和N0分别表示频谱的中心行和中心列。

D0=40; delta=D0;: 设置截止频率距离圆点的距离D0，并将其赋值给delta，delta表示高斯曲线的扩散程度。

for x=1:M: 开始遍历频谱中的每一个像素点的行坐标。

for y=1:N: 在当前行中遍历每一个像素点的列坐标。

%计算点（x,y）到中心点的距离: 计算当前像素点(x,y)到频谱中心点(M0,N0)的距离的平方，并存储在变量d2中。

%计算高斯滤波器: 根据高斯滤波器的定义，使用高斯函数计算当前像素点处的滤波器值，并存储在滤波器矩阵h的对应位置(x,y)中。

end: 结束内层循环。

end: 结束外层循环。

res=h.*Y;: 将计算得到的高斯滤波器h与频谱Y相乘，得到滤波后的频谱结果，并存储在变量res中。

res=real(ifft2(ifftshift(res)));: 对滤波后的频谱进行逆变换，恢复到空域图像，并将结果存储在res中。ifftshift函数将频谱移动回原来的位置，ifft2函数进行二维逆傅里叶变换，real函数取实部。

subplot(223),imshow(res); title('高斯低通滤波所得图像');: 在图像网格中创建第三个子图，显示经过高斯低通滤波后的图像，并设置标题为"高斯低通滤波所得图像"。

subplot(224),imshow(h); title("高斯低通滤波器图象");: 在图像网格中创建第四个子图，显示生成的高斯低通滤波器图像，并设置标题为"高斯低通滤波器图象"。

4.3 算法3——巴特沃斯低通滤波

4.3.1 程序代码

clc
clear
close all
%%
%巴特沃斯低通
figure(3);
I=imread('fingerprint.tif');
subplot(221),imshow(I);
title('原始图像');
Y=fft2(im2double(I));%傅里叶变换
Y=fftshift(Y);%频谱搬移,直流分量搬移到频谱中心
subplot(222), imshow(log(abs(Y)+1),[]); 
title('傅里叶变换取对数所得频谱图像');
[M,N]=size(Y);%获得图像的高度和宽度
h=zeros(M,N);%滤波器函数
%图像中心点
M0=M/2;
N0=N/2;
d0=40;
%巴特沃斯滤波器的阶数
n_0=2;
for x=1:M
    for y=1:N
        distance=sqrt((x-M0)^2+(y-N0)^2);
        h(x,y)=1/(1+(distance/d0)^(2*n_0));
    end
end
%滤波后结果
res=h.*Y;
res=real(ifft2(ifftshift(res)));
subplot(223),imshow(res);
title('巴特沃斯低通滤波所得图像'); 
subplot(224),imshow(h);
title("巴特沃斯低通滤波器图象");