一、实验目的
1.掌握数字图像傅里叶变换；
2.掌握基于频率域的图像滤波算法；
3.能够基于MATLAB实现基于频率域的图像滤波；
4.能够根据实验结果分析各种算法的特点以及应用场合，培养处理实际图像的能力。
二、实验要求
1.实验课前需要写预习实验报告，内容为本次实验要求中的所有程序清单。
2.实验课对预习报告中的编程代码进行上机调试，完成实验指导书中全部实验要求内容。
3.实验课后写出实验报告。报告要求有实验目的，实验内容与步骤，调试完成的准确编程代码，实验小结，回答问题。
三、实验内容（每一个内容编写一个*.m文件）
1.用MATLAB函数对图像(house.tif)完成傅立叶变换；
（1）读入，在图形窗口显示输入图像；
代码：
%显示输入图像
h=imread("house.tif");
figure;
imshow(h);
title("输入图像")结果：
（2）完成傅立叶变换，显示变换后的频谱图；
代码：
%完成傅立叶变换，显示变换后的频谱图
h = im2double(h);
f=fft2(h);
figure;
subplot(1,2,1),imshow(h);
title("原图")
subplot(1,2,2),imshow(uint8(abs(f)));
title("傅里叶变换后的频谱图")结果：
 |
（3）对频谱图采用对数变换，将变换后的图像在图形窗口显示；
代码：
%对频谱图采用对数变换，将变换后的图像在图形窗口显示
fd=log(abs(f)+1); %取模并进行缩放
figure;
subplot(1,2,1),imshow(h);
title("原图")
subplot(1,2,2),imshow(fd,[]);
title("对数变换后的频谱图")结果：
 |
（4）对变换后的频谱进行频移，将零频率部分移到图形的中心位置，并在图形窗口显示；
代码：
%对变换后的频谱进行频移，将零频率部分移到图形的中心位置，并在图形窗口显示
fz=fftshift(fd);
figure;
subplot(1,2,1),imshow(h);
title("原图")
subplot(1,2,2),imshow(abs(fz));
title("频移后的频谱图")结果：
 |
（5）对频谱图采用对数变换，将变换后的图像在图形窗口显示；
代码：
%对频谱图采用对数变换，将变换后的图像在图形窗口显示
fzd=log(abs(fz)+1);
figure;
subplot(1,2,1),imshow(h);
title("原图")
subplot(1,2,2),imshow(fzd,[]);
title("对数变换后的频谱图(fzd)")结果：
 |
（6）获取傅里叶变换后的相位谱图，并在图形窗口显示；
代码：
%获取傅里叶变换后的相位谱图，并在图形窗口显示
xw=angle(f)*180/pi;
figure;
subplot(1,2,1),imshow(h);
title("原图")
subplot(1,2,2),imshow(xw,[]);
title("傅里叶变换后的相位谱图")结果：
 |
（7）对相位谱图采用对数变换，将变换后的图像在图形窗口显示；
代码：
%对相位谱图采用对数变换，将变换后的图像在图形窗口显示
xwd=log(abs(xw)+1);
figure;
subplot(1,2,1),imshow(h);
title("原图")
subplot(1,2,2),imshow(xwd,[]);
title("对数变换后的频谱图")结果：
 |
（8）通过傅里叶反变换得到输入的图像，显示反变换得到的图像。
代码：
%通过傅里叶反变换得到输入的图像，显示反变换得到的图像
fzy=fftshift(f);
fr=real(ifft2(ifftshift(fzy))); %频率域反变换到空间域，并取实部
ff=im2uint8(mat2gray(fr)); %更改图像类型
figure;
subplot(1,2,1),imshow(h);
title("原图")
subplot(1,2,2),imshow(ff);
title("逆傅里叶变换")结果：
 |
2.用MATLAB生成如下的两幅图像（图a：黑白线竖线8条，图b: 黑白线横线32条），图像的大小为
，对这两幅图像进行傅里叶变换，并分别进行频移和对数变换，在一个图形窗口显示：原图像、傅里叶变换后的频谱图和相位图、频移后的频谱图及相位图。
代码：
clc,
clear
lines=8;
m=16; %黑白线宽度
n=256; %图像高度
S=zeros(1,2*m);
S(1,1:m)=ones(1,m);
a=repmat(S,n,lines);
subplot(2,2,1),imshow(a);
title("图a")
lines=32;
m=4; %黑白线宽度
n=256; %图像高度
S=zeros(2*m,1);
S(1:m,1)=ones(1,m);
b=repmat(S,lines,n);
subplot(2,2,2),imshow(b);
title("图b")
%对图a进行傅里叶变换，并分别进行频移和对数变换
%傅里叶变换
fa=fft2(a);
%傅里叶变换的相位图
faxw=angle(fa)*180/pi;
%对数变换
fad=log(abs(fa)+1);
%频移
fp=fftshift(fa);
%频移的相位图
fpxw=angle(fp)*180/pi;
figure;
subplot(2,3,1),imshow(a);
title("图a")
subplot(2,3,2),imshow(uint8(abs(fa)));
title("傅里叶变换后的频谱图")
subplot(2,3,3),imshow(faxw,[]);
title("傅里叶变换的相位图")
subplot(2,3,4),imshow(fad);
title("对数变换后的频谱图")
subplot(2,3,5),imshow(abs(fp));
title("频移后的频谱图")
subplot(2,3,6),imshow(fpxw,[]);
title("频移后的相位图")
%对图b进行傅里叶变换，并分别进行频移和对数变换
%傅里叶变换
fb=fft2(b);
%傅里叶变换的相位图
fbxw=angle(fb)*180/pi;
%对数变换
fbd=log(abs(fb)+1);
%频移
fpb=fftshift(fb);
%频移的相位图
fpbxw=angle(fpb)*180/pi;
figure;
subplot(2,3,1),imshow(b);
title("图b")
subplot(2,3,2),imshow(uint8(abs(fb)));
title("傅里叶变换后的频谱图")
subplot(2,3,3),imshow(fbxw,[]);
title("傅里叶变换的相位图")
subplot(2,3,4),imshow(fbd);
title("对数变换后的频谱图")
subplot(2,3,5),imshow(abs(fpb));
title("频移后的频谱图")
subplot(2,3,6),imshow(fpbxw,[]);
title("频移后的相位图")结果：
生成的图a和图b：
 |
a的原图像、傅里叶变换后的频谱图和相位图、频移后的频谱图及相位图：
 |
b的原图像、傅里叶变换后的频谱图和相位图、频移后的频谱图及相位图：
 |
3.比较傅里叶变换后图像的频谱图和相位图所起的作用。
（1）读入并显示两幅不同的彩色人体图像（man.png和woman.png），大小相同；
（2）分别对上述两幅图像进行傅里叶变换并频移，得到他们的频谱图和相位图；
（3）根据第一个人员图像的频谱图和第二个人员图像的相位图恢复图像；
（4）根据第一个人员图像的相位图和第二个人员图像的频谱图恢复图像；
（5）比较分析图像与原图的差别。
代码：
clc;
clf;
I=imread('man.png');
J=imread('woman.png');
I = rgb2gray(I);
J = rgb2gray(J);%转灰度图像
[M,N]=size(I);
J=imresize(J,[M,N]);%限制图像大小，相位交换需要两个图像大小一致
IF1 = fft2(I);%图像进行傅里叶变换
JF1 = fft2(J);%图像进行傅里叶变换
IF2 = fftshift(abs(IF1));%图像频谱中心化
JF2 = fftshift(abs(JF1));%图像频谱中心化
IABS = abs(IF1);
IAG = angle(IF1);%幅度谱，相位谱
JABS = abs(JF1);
JAG = angle(JF1);
I2=fftshift(abs(IABS)); %log(abs(I2)+1)等于log(abs(IF2)+1)
J2=fftshift(abs(JABS));
figure(1)
subplot(2,4,1);
imshow(I);
title('原始图像1');
subplot(2,4,2);
imshow(log(IABS+1),[]);
title('图像1的幅度谱');
subplot(2,4,3);
imshow(IAG,[]);
title('图像1的相位谱');
subplot(2,4,5);
imshow(J);title('原始图像2');
subplot(2,4,6);
imshow(log(JABS),[]);
title('图像2的幅度谱');
subplot(2,4,7);
imshow(JAG,[]);
title('图像2的相位谱');
subplot(2,4,4);
imshow(log(abs(I2)+1),[]);
title('图像1的频谱中心化');
subplot(2,4,8);
imshow(log(abs(J2)+1),[]);
title('图像2的频谱中心化');
II=ifft2(IF1);
JJ=ifft2(JF1);
figure(2)
subplot(1,2,1);
imshow(abs(II),[]);
title('图像1逆变换');
subplot(1,2,2);
imshow(abs(JJ),[]);
title('图像2逆变换');
%相位交换，注意i，两图像一定要大小一样
ii=IABS.*cos(JAG)+IABS.*sin(JAG).*1i;
jj=JABS.*cos(IAG)+JABS.*sin(IAG).*1i;
iI=abs(ifft2(ii));
jJ=abs(ifft2(jj));
figure(3)
subplot(1,2,1);
imshow(iI,[]);
title('图像1 幅度谱 与 图像2相位谱');
subplot(1,2,2);
imshow(jJ,[]);
title('图像2 幅度谱 与 图像1相位谱')结果：
 |
4.低通滤波和高通滤波的效果比较
（1）读入并显示图像Fig0313(a).tif；
（2）分别对图像进行理想低通滤波和高斯低通滤波，D0=20，60，将以上结果在一个图形窗口显示；
（3）分别对图像进行理想高通滤波和高斯高通滤波，D0=20，60，将以上结果在一个图形窗口显示；
（4）比较分析滤波结果的区别。
区别：
低通滤波（理想/高斯）：
理想低通滤波可能导致明显的振铃效应（Gibbs现象），即在截止频率附近产生明显的震荡和边缘效应。
高斯低通滤波通常具有更平滑的过渡区域，减少了振铃效应的产生。
随着D0值的增加，更多的高频分量被滤除，图像变得更加模糊和平滑。
高通滤波（理想/高斯）：
理想高通滤波可能导致明显的边缘增强和噪声放大效应，突出了图像中的高频细节和噪声。
高斯高通滤波通常具有更平滑的过渡区域，减少了边缘增强和噪声放大的程度。
随着D0值的增加，更少的高频分量被保留，图像中的细节和纹理变得更加模糊和不明显。
代码：
clear;
%读入并显示图像Fig0313(a).tif
f=imread("Fig0313(a).tif");
figure;
imshow(f);
title("原图")
%分别对图像进行理想低通滤波和高斯低通滤波，D0=20，60，将以上结果在一个图形窗口显示
%理想低通滤波,D0=20
lxd=lx_low(f);
%高斯低通滤波,D0=60
gsd=gs_low(f);
figure;
subplot(1,2,1),imshow(lxd);
title("理想低通滤波")
subplot(1,2,2),imshow(gsd);
title("高斯低通滤波")
%分别对图像进行理想高通滤波和高斯高通滤波，D0=20，60，将以上结果在一个图形窗口显示
%理想高通滤波,D0=20
lxg=lx_high(f);
%高斯高通滤波,D0=60
gsg=gs_high(f);
figure;
subplot(1,2,1),imshow(lxg);
title("理想高通滤波")
subplot(1,2,2),imshow(gsg);
title("高斯高通滤波")
function result=lx_low(f)
s=fftshift(fft2(f));
[a,b]=size(s);
a0=round(a/2);
b0=round(b/2);
d=20;
for i=1:a
for j=1:b
distance=sqrt((i-a0)^2+(j-b0)^2);
if distance<=d
h=1;
else
h=0;
end
s(i,j)=h*s(i,j);
end
end
s=uint8(real(ifft2(ifftshift(s))));
result = s;
end
function result = lx_high(f)
s=fftshift(fft2(f));
[a,b]=size(s);
a0=round(a/2);
b0=round(b/2);
d=20;
for i=1:a
for j=1:b
distance=sqrt((i-a0)^2+(j-b0)^2);
if distance>=d
h=1;
else
h=0;
end
s(i,j)=h*s(i,j);
end
end
s=uint8(real(ifft2(ifftshift(s))));
result = s;
end
function result = gs_low(IA)
[f1,f2]=freqspace(size(IA),'meshgrid');
D=60/size(IA,1);
r=f1.^2+f2.^2;
Hd=ones(size(IA));
for i=1:size(IA,1)
for j=1:size(IA,2)
t=r(i,j)/(D*D);
Hd(i,j)=exp(-t);
end
end
Y=fft2(double(IA));
Y=fftshift(Y);
Ya=Y.*Hd;
Ya=ifftshift(Ya);
Ia=real(ifft2(Ya));
result = Ia;
end
function result = gs_high(image)
d0=60; %阈值
[M ,N]=size(image);
img_f = fft2(double(image));%傅里叶变换得到频谱
img_f=fftshift(img_f); %移到中间
m_mid=floor(M/2);%中心点坐标
n_mid=floor(N/2);
h = zeros(M,N);%高斯低通滤波器构造
for i = 1:M
for j = 1:N
d = ((i-m_mid)^2+(j-n_mid)^2);
h(i,j) = 1-exp(-(d)/(2*(d0^2)));
end
end
img_lpf = h.*img_f;
img_lpf=ifftshift(img_lpf); %中心平移回原来状态
img_lpf=uint8(real(ifft2(img_lpf))); %反傅里叶变换,取实数部分
result = img_lpf;
end |
5.（选做）空间域滤波和频率域滤波的比较
（1）读入图像house.tif，对其进行傅里叶变换，显示原图像和傅里叶变换后的频谱图像；
（2）用MATLAB函数fspecial()生成空间滤波器（sobel，laplacian），对图像进行滤波，显示滤波后的图像；
（3）用MATLAB函数将空间滤波器转化为频率域滤波器，对图像进行频率域的滤波；
（4）比较空间滤波和频率域滤波速度和滤波后图像的区别。
区别：
空间域滤波：由于直接在像素上进行操作，可能会对图像的细节造成一定的损失。
频率域滤波：在频率域进行滤波时，可以针对不同的频率成分进行独立的处理，因此能够更好地保护图像的细节和边缘信息。
代码：
clear all;
close all;
figure(1)
h = imread('house.tif');
f = fft2(h);
subplot(2,2,1),imshow(h),title('original image')
subplot(2,2,2),imshow(log(abs(f)+1),[]),title('Fourier Transform1')
% 5-2 调用matlab滤波函数进行滤波增强
h1 = fspecial('laplacian',0);%创建滤波器
h2 = fspecial('sobel');
tic
g1 = imfilter(h,h1,'replicate');%空间域滤波
g2 = imfilter(h,h2,'replicate');
toc %计算空间域滤波时间
disp(['空间域滤波运行时间: ',num2str(toc)]);
subplot(2,2,3),imshow(h-g1),title('laplacian filter函数空间域滤波增强')
subplot(2,2,4),imshow(h-g2),title('sobel filter函数滤波空间域增强')
% 5-3/5-4
figure(2)
subplot(2,3,1);
freqz2(h1);%freqz2()计算线性系统的频率响应,包括幅频响应和相频响应
title('laplacian filter') %查看laplacian滤波器形状--频域
subplot(2,3,4);
freqz2(h2); %查看sobel滤波器形状--频域
title('sobel')
PQ = paddedsize(size(h)); %paddedsize函数为自定义函数，计算填充尺寸以供基于FFT的滤波器
H1 = freqz2(h1,PQ(1),PQ(2));
H2 = freqz2(h2,PQ(1),PQ(2));
H11 = ifftshift(abs(H1)); %频移，得到频域滤波函数频谱。
H22 = ifftshift(abs(H2));
subplot(2,3,2),imshow(abs(H1),[]) %显示频谱
subplot(2,3,3),imshow(abs(H11),[]) %频移后，显示频谱
subplot(2,3,5),imshow(abs(H2),[])
subplot(2,3,6),imshow(abs(H22),[])
figure(3)
tic
gf1 = dftfilt(h, H11); %dftfilt函数为自定义函数，实现傅里叶变换并滤波
gf2 = dftfilt(h, H22); %频域滤波
toc %计算频域滤波时间
disp(['频域滤波运行时间: ',num2str(toc)]);
subplot(2,3,[1,4]),imshow(h),title('原图')
subplot(2,3,2),imshow(g1,[]),title('空间域滤波（laplacian）')
subplot(2,3,3),imshow(g2,[]),title('空间域滤波（sobel）')
subplot(2,3,5),imshow(gf1,[]),title('频率域滤波（laplacian）')
subplot(2,3,6),imshow(gf2,[]),title('频率域滤波（sobel）')
figure(4)
subplot(2,3,[1,4]),imshow(h),title('原图')
subplot(2,3,2),imshow(h-g1,[]),title('空间域滤波（laplacian）增强')
subplot(2,3,3),imshow(h-g2,[]),title('空间域滤波（sobel）增强')
subplot(2,3,5),imshow(h-uint8(gf1),[]),title('频率域滤波（laplacian）增强')
subplot(2,3,6),imshow(h-uint8(gf2),[]),title('频率域滤波（sobel）增强')
function PQ = paddedsize(AB, CD)
% 计算填充尺寸以供基于FFT的滤波器
% PQ = PADDEDSIZE(AB),AB = [A B], PQ = 2 * AB
%
% PQ = PADDEDSIZE(AB, 'PWR2')， PQ（1） = PQ（2） = 2 ^ nextpow2(2 * m), m =
% MAX(AB).
%
% PQ = PADDEDSIZE(AB, CD)，AB = [A B], CD = [C D]
%
% PQ = PADDEDSIZE(AB, CD, 'PWR2')， PQ（1） = PQ（2） = 2 ^ nextpow2(2 * m), m =
% MAX([AB CD]).
if nargin == 1
PQ = 2 * AB;
elseif nargin == 2 & ~ischar(CD)
PQ = AB + CD -1;
PQ = 2 * ceil(PQ / 2); % ceil（N）返回比N大的最小整数，为了避免出现奇数，因为处理偶数数组快
elseif nargin == 2
m = max(AB);
P = 2 ^ nextpow2(2 * m); % nextpow2（N）返回第一个P，使得2. ^ P> = abs（N）。
% 对于FFT运算，找到最接近两个序列长度的2 的幂次方通常很有用。
PQ = [P, P];
elseif nargin == 3
m = max([AB CD]);
P = 2 ^ nextpow2(2 * m);
PQ = [P, P];
else
error('Wrong number of input')
end
end
function g = dftfilt(f, H)
% DFTFILT performs frequency domain filtering.
% G = DFTFILT(F, H) filters F in the frequency domain using the filter
% transfer function H. The output, G, is the filtered image, which has
% the same size as F. DFTFILT automatically pads F to be the same size as
% H. Function PADEDESIZE can be used to determine an appropriate size
% for H.
% G = DFTFILT（F，H）使用滤波器传递函数H在频域中对输入图像F滤波。 输出G是滤波后的图像，
% 其大小与F相同。DFTFILT自动将F填充到与H相同的大小 ，PADEDESIZE函数可用于确定H的合适大小。
%
% DFTFILT assumes that F is real and that H is a real, uncentered,
% circularly-symmetric filter function.
% DFTFILT假设F是实数，H是实数，未中心，循环对称的滤波函数。
% Obtain the FFT of the padded input.
% 获取填充之后的FFT变换
F = fft2(f, size(H, 1), size(H, 2));
% Perform filtering
% 滤波
g = real(ifft2(H .* F));
% Crop to orihinal size
% 剪切到原始尺寸
g = g(1 : size(f, 1), 1 : size(f, 2));
end结果：
 |
四、问题与讨论
1.为什么要在频率域对图像进行滤波？
在频率域对图像进行滤波的原因有以下几点：
可以利用频率成分和图像外表之间的对应关系。一些在空间域表述困难的增强任务，在频率域中变得非常普通。
滤波在频率域更为直观，它可以解释空间域滤波的某些性质。
可以在频率域指定滤波器，做反变换，然后在空间域使用结果滤波器作为空间域滤波器的指导。
- 根据题1的实验结果分析（1）移频前后的频谱图有什么区别？（2）对数变换对频谱图的影响怎样？
根据题1的实验结果：
（1）移频前后的频谱图的区别在于频率轴上的位置发生了变化。具体来说，移频后的频谱图整体向右移动，而移频前的频谱图则整体向左移动。这是因为移频操作实质上是在时域内将信号的相位滞后或提前，这在频谱图上就表现为频率的左移或右移。
（2）对数变换对频谱图的影响主要体现在将频谱图压缩和拉伸，但不会改变各频率成分的相对大小。具体来说，对数变换可以使得低频成分在频谱图上得到增强，而高频成分则相对减弱。这有利于我们在处理图像等信号时，更清楚地观察到低频成分的细节和特征。
3.根据题3的实验结果，分析频谱图和相位图的作用。
根据题3的实验结果，频谱图和相位图在信号处理和分析中起着至关重要的作用。
频谱图是一个表示信号频率内容及其强度的图表。在频谱图中，每个频率分量的强度用不同颜色的条纹表示。通过观察频谱图，我们可以了解信号中各种频率分量的分布情况，以及各个频率分量之间的相对强度。这对于分析信号的频率特性、噪声情况以及进行频率域滤波等处理非常有帮助。
相位图则是一个表示信号相位关系的图表。在相位图中，每个点的颜色和位置表示信号在特定时刻的相位状态。通过观察相位图，我们可以了解信号的相位变化情况，例如是否存在相位跳跃、相位偏移等现象。这对于分析信号的稳定性、预测信号的变化趋势以及进行相位域处理等非常有帮助。
4.根据题4的实验结果，分别分析（1）低通滤波和高通滤波对图像的影响，（2）截止频率对低通滤波的影响，（3）相同截止频率时，理想低通滤波器和高斯低通滤波器对图像的影响有什么不同。
根据题4的实验结果可得
（1）低通滤波和高通滤波对图像的影响：
低通滤波器通常用于保留图像的平滑部分，去除高频噪声，例如边缘和细节。这可以使得图像看起来更加模糊，因为高频成分被抑制。高通滤波器则相反，它能够强化图像的边缘和细节，去除平滑区域中的噪声。因此，高通滤波后的图像会显得更加清晰，突出了高频成分。
（2）截止频率对低通滤波的影响：
低通滤波器的性能受到其截止频率的影响。当提高截止频率时，更多的高频成分会被保留，导致图像的细节和边缘更加明显。这可能会导致图像看起来更加清晰，但同时也可能增强图像中的噪声。相反，当降低截止频率时，更多的高频噪声会被抑制，导致图像更加平滑，但细节和边缘可能会变得模糊。因此，选择合适的截止频率对于实现最佳的低通滤波效果非常重要。
（3）相同截止频率时，理想低通滤波器和高斯低通滤波器对图像的影响有什么不同：
理想低通滤波器在频率域中完全抑制了高于截止频率的频率分量，而在其他区域则完全通过。这种类型的滤波器对于去除高频噪声非常有效，但可能会导致边缘的锐化程度不够。高斯低通滤波器则采用平滑过渡的方式抑制高频成分，而不是完全抑制。这可以更好地保护图像的边缘，避免出现明显的边缘效应，但同时也可能会导致平滑区域中的一些细节丢失。因此，在选择使用理想低通滤波器还是高斯低通滤波器时，需要根据具体的应用场景和需求进行权衡。
5.（选做） 根据题5的实验结果，比较空间滤波和频率域滤波对图像的影响。
空间滤波和频率域滤波对图像的影响主要有以下区别：
空间滤波是在图像的像素域上直接对图像进行处理，通过在像素级别上应用某些操作，以达到增强或改进图像质量的目的。空间滤波器通常由线性或非线性滤波器组成，线性滤波器包括均值滤波、高斯滤波等，非线性滤波器包括中值滤波、边缘保持滤波等。不同的空间滤波器具有不同的特性，可以对图像进行平滑、锐化、边缘检测等处理。
频率域滤波则是在图像的频率域上对图像进行处理，通过改变图像的频率成分来达到增强或改进图像质量的目的。频率域滤波器通常包括低通滤波器、高通滤波器、带通滤波器等。低通滤波器可以去除图像中的高频噪声，保留图像的平滑部分；高通滤波器则可以强化图像的边缘和细节，去除平滑区域中的噪声；带通滤波器则可以根据需要选择性地保留图像中的某些频率成分。
扫一扫
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
