图像处理学习笔记08：频域滤波

计算与观察二维DET

在MATLAB中傅里叶变换可以是使用快速傅里叶变换（FFT）实现。使用函数fft2就可以实现，语法形式F=fft2(f),这个函数返回的傅里叶变换大小仍为MxN。傅里叶谱可以使用abs函数计算每个元素实部与虚部平方和的平方根获得。
例

f= imread("Fig0303(a).tif");
figure;
subplot(2,3,1);
imshow(f);
%计算f的傅里叶变换并绘制频谱
F=fft2(f);
S=abs(F);
subplot(2,3,2);
imshow(S,[]);
%中心变换，把变换的原点移到频率矩阵中心，fftshift是通过变换f的象限进行操作的。操作结果与将输入图像乘负一的x+y次效果相同
Fc=fftshift(F);
subplot(2,3,3);
imshow(abs(Fc),[]);
%使用灰度变换中的对数变换，将图像中暗的区域扩展
S2=log(1+abs(Fc));
subplot(2,3,4);
imshow(S2,[]);
%计算相角，使用imag和real分别取出矩阵的实部和虚部，使用atan2函数就可以算出相角。也可以直接使用angle函数得到矩阵的相角。
phi=atan2(imag(F),real(F));%phi=angle(F);
subplot(2,3,5);
imshow(phi);
f1=S.*exp(1i*phi);
subplot(2,3,6);
imshow(f1);

DFT滤波的基本步骤

1. 用函数tofloat把输入图像变换为浮点图像。[f,revertclass]=tofloat(f);
2. 用函数paddedsize获得填充参数：PQ=paddedsize(size(f));
3. 得到有填充的傅里叶变换:F=fft2(f,PQ(1),PQ(2));
4. 生成大小为PQ(1)xPQ(2)的滤波函数H，滤波器函数必须是四角凸起的，如果突起居中，令H=fftshift(H)
5. 用滤波器乘以FFT变换G=H.*F;
6. 获得G的逆FFT变换g=ifft2(G);
7. 修剪左上部矩形为原始大小g=g(1:size(f,1),1:size(f,2));
8. 把滤波过的图像变换为输入图像的类g=revertclass(g);

从空域滤波器到频域滤波器

函数freqz2计算FIR滤波器的频率响应，语法形式H=freqz2(h,R,C);h是二维空域滤波器，R是行数，C是我们希望的滤波器H的列数，H是相应的二维频域滤波器。
例 空域和频域滤波器的比较。

f=imread("Fig0309(a).tif");
f=tofloat(f);%将图像转化为浮点类型
figure;
subplot(5,2,1);
imshow(f);
F=fft2(f);%对f进行傅里叶变换
S=fftshift(log(1+abs(F)));%先对F进行灰度变换，然后进行中心变换。
subplot(5,2,2);
imshow(S,[]);
h=fspecial('sobel');
h=h.';
% freqz2(h);
PQ=paddedsize(size(f));
H=freqz2(h,PQ(1),PQ(2));
H1=ifftshift(H);%将中心变换恢复，使原点在频域的左上角
%计算图像幅值，绘制图像
subplot(5,2,3)
imshow(abs(H),[]);
subplot(5,2,4)
imshow(abs(H1),[])
gs=imfilter(f,h);%采用默认值0填充边缘
gf=dftfilt(f,H1);%把滤波后的图像变换为与输入相同的类
%图像中的灰色是由于含有负数引起的
subplot(5,2,5)
imshow(gs,[]);
subplot(5,2,6)
imshow(gf,[]);
%显示出计算图像的绝对值
subplot(5,2,7)
imshow(abs(gs),[]);
subplot(5,2,8)
imshow(abs(gf),[]);
%使用阈值建立一幅二值图像，以便边缘看的更加清晰
subplot(5,2,9)
imshow(abs(gs)>0.2*abs(max(gs(:))));
subplot(5,2,10)
imshow(abs(gf)>0.2*abs(max(gf(:))));
% figure;
% h=fftshift(h);
% freqz2(h);

频域低通（平滑）滤波器

clear;clc;
f=imread("Fig0313(a).tif");
[f,revertclass]=tofloat(f);
PQ=paddedsize(size(f));%填充参数
[U,V]=dftuv(PQ(1),PQ(2));%距离计算
D=hypot(U,V);%求平方和再开方
D0=0.05*PQ(2);%截止频率
F=fft2(f,PQ(1),PQ(2));%滤波后的傅里叶谱
H=exp(-(D.^2)/(2*(D0^2)));%高斯低通滤波器
g=dftfilt(f,H);%对图像进行高斯低通滤波处理
g=revertclass(g);
figure;
subplot(2,2,1);
imshow(f);
subplot(2,2,2);
imshow(fftshift(H));
subplot(2,2,3);
imshow(log(1+abs(fftshift(F))),[]);
subplot(2,2,4);
imshow(g);

高通（锐化）频域滤波器

例 绘制三种高通滤波器的图像

%分别为理想高通滤波器、巴特沃斯滤波器、高斯滤波器
H=fftshift(hpfilter('ideal',500,500,50));
H1=fftshift(hpfilter('btw',500,500,50));
H2=fftshift(hpfilter('gaussian',500,500,50));
mesh=(double(H(1:10:500,1:10:500)));
axis tight;
colormap([0 0 0]);
axis off;
figure;
subplot(1,3,1);
imshow(H,[]);
subplot(1,3,2);
imshow(H1,[]);
subplot(1,3,3);
imshow(H2,[]);

例 高通滤波处理图像

f=imread("Fig0318(a).tif");
figure;
subplot(1,2,1);
imshow(f);
PQ=paddedsize(size(f));%获得填充参数
D0=0.05*PQ(1);%截止频率
H=hpfilter('gaussian',PQ(1),PQ(2),D0);
g=dftfilt(f,H);%使用高斯高通滤波器对图像进行处理
subplot(1,2,2);
imshow(g);

图像中的边缘和其他锐利转换的灰度增强了。

高频强调滤波

高通滤波器偏移了dc的0项，因此减少了图像中平均值为0的值。一种补偿方法就是为高通滤波器加上偏移量。如果偏移量与将滤波器乘以某个大于一的常数结合起来，因为这个乘数突出了高频部分，这种方法就叫高频强调滤波。偏移量要小于乘数，低频增强的影响就弱于高频增强的影响。
例 高频强调滤波和直方图均衡结合

f=imread("Fig0319(a).tif");
PQ=paddedsize(size(f));
D0=0.05*PQ(1);
HBW=hpfilter('btw',PQ(1),PQ(2),D0,2);
H=0.5+2*HBW;
gbw=dftfilt(f,HBW,'fltpoint');
gbw=gscale(gbw);
ghf=dftfilt(f,H,'fltpoint');
ghf=gscale(ghf);
ghe=histeq(ghf,256);
figure;
subplot(2,2,1);
imshow(f);
subplot(2,2,2);
imshow(gbw);
subplot(2,2,3);
imshow(ghf);
subplot(2,2,4);
imshow(ghe);

选择性滤波

例 用陷波滤波器减少波纹模式

clear;
clc;
f=imread('Fig0321(a).tif');
figure;
subplot(2,3,1);
imshow(f);
[M,N]=size(f);
[f,revertclass]=tofloat(f);
F=fft2(f);
S=gscale(log(1+abs(fftshift(F))));
subplot(2,3,2);
imshow(S);
%使用imtool函数交互式地确定峰值坐标
C1=[99 154;128 163];
%陷波滤波器
H1=cnotch('gaussian','reject',M,N,C1,5);
%计算滤波后的频谱，并显示图像
P1=gscale(fftshift(H1).*(tofloat(S)));
subplot(2,3,3);
imshow(P1);
g1=dftfilt(f,H1);
g1=revertclass(g1);
subplot(2,3,4);
imshow(g1);
%重复下列c2，以减小高频干扰成分
C2=[99 154;128 163;49 160; 133 233; 55 132;108 225;112 74];
H2=cnotch('gaussian','reject',M,N,C2,5);
%计算滤波后的频谱，并显示图像
P2=gscale(fftshift(H2).*(tofloat(S)));
subplot(2,3,5);
imshow(P2);
g2=dftfilt(f,H2);
g2=revertclass(g2);
subplot(2,3,6);
imshow(g2);

例 用陷波滤波器减少因有故障的成像设备而导致的周期干扰。

clear;clc;
f=imread("Fig0322(a).tif");
imshow(f);
[M,N]=size(f);
[f,revertclass]=tofloat(f);
F=fft2(f);
S=gscale(log(1+abs(fftshift(F))));
imshow(S);%垂直方向上的能量脉冲是由干扰造成的
H=recnotch('reject','vertical',M,N,3,15,15);
imshow(fftshift(H));
g=dftfilt(f,H);
g=revertclass(g);
imshow(g);
%用陷波滤波器代替带阻滤波器，在垂直轴上隔离干扰的频率。
Hrecpass=recnotch('pass','vertical',M,N,3,15,15);
interference=dftfilt(f,Hrecpass);
imshow(fftshift(Hrecpass));
interference=gscale(interference);
imshow(interference);