Matlab Chap3 频率域滤波

主要的reference：《数字图像处理 MATLAB 冈萨雷斯》

1. 二维傅里叶变换

注：常用数学符号的 LaTeX 表示方法

F(u,v)=∑M−1x=0∑N−1y=0f(x,y)e−j2π(ux/M+yv/N) (DFT)

f(x,y)=1MN∑M−1u=0∑N−1v=0F(u,v)ej2π(ux/M+vy/N) (IDFT)

F(0,0)为直流分量，F(0,0)=1MN∑M−1u=0∑N−1v=0f(x,y) F(0,0)是f(x,y)平均值的MN倍。

功率谱P(u,v)=F(u,v)2 .

1.若f(x,y)是实函数，F(u,v)仍然是复数，但关于原点共轭对称。F(u,v) = F*(-u,-v);

2.频谱（幅度）也关于原点对称F(u,v) = |F(-u,-v)|

3.F(u,v) = F(u+k1M,v+k2N)。DFT在u，v方向上每隔M,N个长度周期重复。

同理，f(x,y)（IDFT）也是每隔M，N周期重复。

4.fftshift=DFT−1x+yf(x,y)=F(u−M/2,v−N/2) 频移特性，可使得频谱中心在矩阵中心，而非第一个元素。

5.频谱拉伸：举例，lena.bmp图片的第一个通道G分量（0~255）其频谱fft2(f)=F的范围为（1.9034,~48159521）。对于imshow(abs(F),[])是把F的范围扩展到0`255，显然截断了大部分信息。故使用log(1+abs(F))变为（0~17.69），再用imshow(,[])灰度拉伸到0~255.

2. 填零补充问题

由于离散傅里叶变换是，周期。空域卷积是循环卷积，在做DFT的时候，会导致上边沿模糊，而左右边沿不模糊。是因为没有填零周期复制。为了填0，采用填充处理。如书中图所示，会导致左右边沿，上边沿均模糊。而与此类似的处理，是直接imfilter(f,h)，直接在空域做。频域做：如下代码。可以直接把如下代码集成为一个子函数，见：

g = dftfil(f,H,classout);%classout = 'original'or'fltpoit'

f = imread('lena.bmp');f = f(:,:,1);F = fft2(f);PQ=paddedsize(size(f));Fp = fft2(f,PQ(1),PQ(2));Hp = lpfilter('gaussian',PQ(1),PQ(2),2*sig); % PQ=[512 512]% figure;imshow(fftshift(Hp),[])Gp = Hp.*Fp;gp = ifft2(Gp);gpc = gp(1:size(f,1), 1:size(f,2));%截取

chap3.4从空间滤波器到频域滤波器

关键代码：freqz(h);freqz(h,PQ(size(f,1), PQ(size(f,2))));dftfilt(f,H1);

f = imread('lena.bmp');f = f(:,:,1);h = fspecial('sobel')'; % 增强垂直边缘%      1     0    -1%      2     0    -2%      1     0    -1figure(1),freqz2(h); title('与空域Sobel对应的频率中的H');% uses [n2 n1] = [64 64].% size_h = size(temp)PQ = paddedsize(size(f));H = freqz2(h,PQ(1),PQ(2)); % 产生的滤波器原点在矩阵中心处H1 = ifftshift(H); % 迁移原点到左上角figure,mesh(abs(H1(1:10:512,1:10:512)));%生成滤波后的图像，空域与频域处理相比较gs = imfilter(double(f),h);gf = dftfilt(f,H1);%零填充处理figure(1);imshow(gs,[]);%gs,gf中有大量负值figure(2);imshow(gf,[]);figure(3);imshow(abs(gs),[]);figure(4);imshow(abs(gf),[]);figure(5);imshow(abs(gs) > 0.2*abs(max(gs(:))))%阈值化了，只显示比阈值大的像素figure(6);imshow(abs(gf) > 0.2*abs(max(gf(:))))d = abs(gs-gf);max(d(:))min(d(:))%由d可知，频域处理结果与空域处理结果，是一样的。

chap3.5在频率域中直接生成滤波器

距离函数

实现循环对称滤波器，首先是距离函数。使得到原点的距离相等的点具有相同的距离。[U,V] = dftuv(M,N);就是实现这样的功能。只是原点在左上角。

然后用距离函数生成别的滤波器函数。

使用lpfilter生成滤波器函数(频率域的中心在左上角的) H = lpfilter(type,M,N,D0,n) 把函数集成为了而已，没有什么高深的代码。

type =‘ideal’理想低通滤波器，or ‘btw’巴特沃斯低通滤波器 or ‘gaussian’高斯低通滤波器

D0是截止频率在距离中心为D0的地方；n是巴特沃斯滤波器的阶数，其余两个不需要此参数。

M,N为矩阵大小

F1 = fft2(f); % 注意 F1 和 下面 F 的区别figure(1);imshow(log(1+abs(fftshift(F1))),[])title('傅立叶频谱图像')PQ = paddedsize(size(f));[U V] = dftuv(PQ(1),PQ(2));%生成距离函数D0 = 0.05*PQ(2);%截止频率所在的距离处F = fft2(f,PQ(1),PQ(2));%零填充之后的频谱figure(2);imshow(log(1+abs(fftshift(F))),[]);title('傅立叶频谱图像 零填充')H = exp(-(U.^2+V.^2)/(2*(D0^2)));figure(3);imshow(fftshift(H),[]);title('高斯低通滤波器频谱图像');g = dftfilt(f,H);figure(4);imshow(g,[])title('高斯低通处理后图像')% figure;mesh(fftshift(abs(H(1:20:1024,1:20:1024))));%高斯低筒滤波器的三维网状图

mesh的使用(surf与此类似)

mesh(fftshift(abs(H(1:20:1024,1:20:1024))));%%控制整个图的显示axis([0 D0 0 D0 0 1]);colormap([0 0 0])%颜色变为黑白axis off; grid off;%关闭网格线和坐标轴view(-25,30);view(-25,0);%改变视角（方位角，仰角）

高通滤波的情况

由于hpfilter = 1-lpfilter；故这样得到的滤波器，

最后，所有源代码

f = imread('Fig0405(a)(square_original).tif');F =fft2(f);Fc = fftshift(F);S = log(1+abs(Fc));[M,N]=size(f);%高斯低通滤波器sig = 10;H = lpfilter('gaussian',M,N,sig); % max(H(:)) = 1% imshow(1-H,[]) % 显示表明滤波器图像未置中title('滤波器频谱（取反）图像');G = H.*F;%注：此处g含有负值g = ifft2(G);%G是实对称的，ifft就会检验到G是实对称，则g是实数% imshow(g,[]);title('不使用填充的频域低通滤波处理后的图像')%PQ = 2*size(f)PQ = paddedsize(size(f)); % size(f)=[256 256]%对原始图像的周围进行0填充Fp = fft2(f, PQ(1),PQ(2));Hp = lpfilter('gaussian',PQ(1),PQ(2),2*sig); % PQ=[512 512]% figure;imshow(fftshift(Hp),[])Gp = Hp.*Fp;gp = ifft2(Gp);% figure;imshow(gp,[]);title('使用填充的频域低通滤波处理后的图像')gpc = gp(1:size(f,1), 1:size(f,2));%截取% imshow(gpc,[]);h = fspecial('gaussian',15,7);gs = imfilter(f,h);% imshow(gs,[]);title('使用空间滤波器处理后的图像');% figure,mesh(abs(Hp(1:10:512,1:10:512)));g = dftfilt(f,H);  % 几步合并为一步 % 要求 H 原点在左上角% figure;imshow(g,[]);%%chap3.4从空间滤波器到频域滤波器f = imread('lena.bmp');f = f(:,:,1);h = fspecial('sobel')'; % 增强垂直边缘%      1     0    -1%      2     0    -2%      1     0    -1figure(1),freqz2(h); title('与空域Sobel对应的频率中的H');% uses [n2 n1] = [64 64].% size_h = size(temp)PQ = paddedsize(size(f));H = freqz2(h,PQ(1),PQ(2)); % 产生的滤波器原点在矩阵中心处H1 = ifftshift(H); % 迁移原点到左上角figure,mesh(abs(H1(1:10:512,1:10:512)));%生成滤波后的图像，空域与频域处理相比较gs = imfilter(double(f),h);gf = dftfilt(f,H1);%零填充处理figure(1);imshow(gs,[]);%gs,gf中有大量负值figure(2);imshow(gf,[]);figure(3);imshow(abs(gs),[]);figure(4);imshow(abs(gf),[]);figure(5);imshow(abs(gs) > 0.2*abs(max(gs(:))))%阈值化了，只显示比阈值大的像素figure(6);imshow(abs(gf) > 0.2*abs(max(gf(:))))d = abs(gs-gf);max(d(:))min(d(:))%由d可知，频域处理结果与空域处理结果，是一样的。%%chap3.5在频率域中直接生成滤波器% [U,V] = dftuv(7,5);% % [U,V] = dftuv(8,5);% D = U.^2 + V.^2;% fftshift(D);F1 = fft2(f); % 注意 F1 和 下面 F 的区别figure(1);imshow(log(1+abs(fftshift(F1))),[])title('傅立叶频谱图像')%低通滤波，高斯低通滤波PQ = paddedsize(size(f));[U V] = dftuv(PQ(1),PQ(2));%生成距离函数D0 = 0.05*PQ(2);%截止频率所在的距离处F = fft2(f,PQ(1),PQ(2));%零填充之后的频谱figure(2);imshow(log(1+abs(fftshift(F))),[]);title('傅立叶频谱图像 零填充')H = exp(-(U.^2+V.^2)/(2*(D0^2)));figure(3);imshow(fftshift(H),[]);title('高斯低通滤波器频谱图像');g = dftfilt(f,H);figure(4);imshow(g,[])title('高斯低通处理后图像')% figure;mesh(fftshift(abs(H(1:20:1024,1:20:1024))));%高斯低筒滤波器的三维网状图%%高通滤波的情况PQ = paddedsize(size(f));D0 = 0.05*PQ(1);HBW = hpfilter('btw',PQ(1),PQ(2),D0,2);% mesh(fftshift(HBW(1:40:1024,1:40:1024)));H = 0.5+2*HBW;%偏离了直流项，高频强调滤波gbw = dftfilt(f,HBW);%频域处理% 使用了 gscale(gbw) 之后，imshowMy(gbw) 等价于 imshowMy(gbw,[])gbw = gscale(gbw); figure(1);imshow(gbw,[])title('高通滤波后的图像')gbe = histeq(gbw,256);%直方图均衡化figure(2);imshow(gbe,[])title('高通滤波并经过直方图均衡化后的图像')ghf = dftfilt(f,H);%频域滤波处理ghf = gscale(ghf);figure(3);imshow(ghf,[])title('高频强调滤波后的图像')ghe = histeq(ghf,256);figure(4);imshow(ghe,[])title('高频强调滤波并经过直方图均衡化后的图像')