数字图像处理学习笔记5：频率域滤波1（傅里叶频谱图，低通滤波-平滑，高通滤波-锐化）

---

前言

频率域滤波步骤：

频率域滤波分为低通滤波（平滑）与高通滤波（锐化）：
$$低通滤波=\{\begin{array}{l}理想低通滤波\\ 布特沃斯低通滤波\\ 高斯低通滤波\end{array}$$
$$高通滤波=\{\begin{array}{l}理想高通滤波\\ 布特沃斯高通滤波\\ 高斯高通滤波\end{array}$$

---

一、傅里叶变换：傅里叶频谱图

使用下列代码得到图像的频谱图

I=imread('1.jpg');
I=rgb2gray(I);
figure
imshow(I)
I=im2double(I);
F=fft2(I);
F=fftshift(F);
F=abs(F);
T=log(F+1);
figure
imshow(T,[]);

结果：
若仅使用频谱图来进行简单滤波，只需要了解频谱图中心为0频分量，越接近边缘则是高频分量，其中高频分量代表图像中灰度变化较大的信息，如噪声和细节，低频则相反。

---

比如对上图加入椒盐噪声，得到：
与上面的频谱图比较，发现频谱图外层的亮点明显增多，因为增加了椒盐噪声，高频分量变多。

使用频谱图进行滤波，就是在频谱图中减去低频或高频分量。

---

二、低通滤波

$$低通滤波=\{\begin{array}{l}理想低通滤波\\ 布特沃斯低通滤波\\ 高斯低通滤波\end{array}$$

---

1.理想低通滤波

理想低通滤波是在频谱图中，以0频分量（中心点）为圆心，以r为半径作的圆中，使圆外的所有高频分量置零，即只保留圆内的低频分量，以此去除噪声。

理想低通滤波器函数为：
$$H(u,v)=\{\begin{array}{ll}1,& \text{D<D0}\\ 0,& \text{D>D0}\end{array}$$
式中D为点到中心点的距离,$D_0$为设置半径

函数各种表现形式如下：

---

以上面加了椒盐噪声的图为例进行理想低通滤波：

代码如下：

im1=imread('椒盐噪声图像.jpg');
im1=rgb2gray(im1);
figure
imshow(im1)

r=80;   %设置滤波圆半径参数为80
img_f=fftshift(fft2(double(im1)));  %傅里叶变换得到频谱
[m,n]=size(img_f);
P_x=fix(m/2);  
P_y=fix(n/2);  %获取圆心坐标
H=zeros(m,n);  
for j=1:n
    for i=1:m
        d=sqrt((i-P_x)^2+(j-P_y)^2);   %计算两点之间的距离，判断在圆外还是圆内
        if d<=r
            H(i,j)=1;   %圆内数值不变
        else
            H(i,j)=0;   %圆外置零，即滤除高频分量
        end
    end
end
img_lpf=H.*img_f; 
img_lpf=ifftshift(img_lpf);            %傅里叶反变换
img_lpf=uint8(real(ifft2(img_lpf)));   %取实数部分
figure
imshow(img_lpf)
imwrite(img_lpf,'理想低通滤波.jpg')

---

滤波前后图像结果如下：

---

2.布特沃斯低通滤波

与理想低通滤波的“极端”不同，布特沃斯低通滤波器以渐变的方式改变频率的数值大小。

布特沃斯低通滤波函数为：
$$H(u,v)=\frac{1}{1+[D(u,v)/D_0{]}^{2n}}$$
式中，D(u,v)为点到中心点的距离，$D_0$为设置的半径

各种表现形式为：

---

代码如下：

im1=imread('椒盐噪声图像.jpg');
im1=rgb2gray(im1);
figure
imshow(im1)

r=80;   %设置滤波圆半径参数为80，可调
img_f=fftshift(fft2(double(im1)));  %傅里叶变换得到频谱
[m,n]=size(img_f);
O_x=fix(m/2);  
O_y=fix(n/2);  %获取圆心坐标
img=zeros(m,n);  %提前定义滤波后的频谱，提高运行速度
for j=1:n
    for i=1:m
        d=sqrt((i-O_x)^2+(j-O_y)^2);    %计算两点之间的距离
        H(i,j)=1/((1+d/r)^4);   %布特沃斯滤波器，这里设n=2，即2阶滤波器，n可调

    end
end
img=H.*img_f; %滤波
img=ifftshift(img);            %傅里叶反变换
img=uint8(real(ifft2(img)));   %取实数部分
figure
imshow(img)

---

处理前后图像：

---

3.高斯低通滤波

高斯滤低通波器函数为：
$$H(u,v)=e^{-D^2(u,v)/2D_0^2}$$

各表现形式如下：

---

代码如下：

im1=imread('椒盐噪声图像.jpg');
im1=rgb2gray(im1);
figure
imshow(im1)

r=80;   %设置滤波圆半径参数为80，可调
img_f=fftshift(fft2(double(im1)));  %傅里叶变换得到频谱
[m,n]=size(img_f);
O_x=fix(m/2);  
O_y=fix(n/2);  %获取圆心坐标
img=zeros(m,n);  %提前定义滤波后的频谱，提高运行速度
for j=1:n
    for i=1:m
        d=sqrt((i-O_x)^2+(j-O_y)^2);    %计算两点之间的距离
         H(i,j)=exp((-d.^2)/(2*(r).^2));   %高斯滤波器
        end
end
img=H.*img_f; %滤波
img=ifftshift(img);            %傅里叶反变换
img=uint8(real(ifft2(img)));   %取实数部分
figure
imshow(img)

---

滤波前后结果：

---

4.小结

可以发现，高斯滤波和理想低通滤波程序中，只能调节半径一个参数，属于两个极端，而布特沃斯滤波器不仅可以设置半径，还可以设置阶数，通过调节阶数可以接近理想低通滤波与高斯低通滤波结果，处于两者之间的过渡。

---

三、高通滤波

低通滤波是减少高频分量，同理，高通滤波则是减少低频分量，并且高通滤波器可以由对应低通滤波器给出：
$$H_{HP}(u,v)=1-H_{LP}(u,v)$$

则对应有三个高通滤波器：
$$高通滤波=\{\begin{array}{l}理想高通滤波\\ 布特沃斯高通滤波\\ 高斯高通滤波\end{array}$$

---

1.理想高通滤波

理想高通滤波函数：
$$H(u,v)=\{\begin{array}{ll}0,& \text{D<D0}\\ 1,& \text{D>D0}\end{array}$$
频谱图中心的值与原图像的灰度平均值相关，高通滤波中，将中心点置零，则滤波后的图像平均值变为0，图像表现较暗。

---

代码如下：

im1=imread('椒盐噪声图像.jpg');
im1=rgb2gray(im1);
figure
imshow(im1)

r=80;   %设置滤波圆半径参数为80
img_f=fftshift(fft2(double(im1)));  %傅里叶变换得到频谱
[m,n]=size(img_f);
P_x=fix(m/2);  
P_y=fix(n/2);  %获取圆心坐标
img_lpf=zeros(m,n);  %提前定义滤波后的频谱，提高运行速度
for j=1:n
    for i=1:m
        d=sqrt((i-P_x)^2+(j-P_y)^2);   %计算两点之间的距离，判断在圆外还是圆内
        if d<=r
            H(i,j)=0;   %圆内置零，去除低频分量
        else
            H(i,j)=1;   %圆外不变
        end
        img_lpf(i,j)=H(i,j)*img_f(i,j);
    end
end
 
img_lpf=ifftshift(img_lpf);            %傅里叶反变换
img_lpf=uint8(real(ifft2(img_lpf)));   %取实数部分
figure
imshow(img_lpf)
imwrite(img_lpf,'理想高通滤波.jpg')

---

2.布特沃斯高通滤波

布特沃斯滤波函数为：
$$H(u,v)=\frac{1}{1+[D_0/D(u,v){]}^{2n}}$$

---

代码如下：

im1=imread('椒盐噪声图像.jpg');
im1=rgb2gray(im1);
figure
imshow(im1)

r=80;   %设置滤波圆半径参数为80，可调
img_f=fftshift(fft2(double(im1)));  %傅里叶变换得到频谱
[m,n]=size(img_f);
O_x=fix(m/2);  
O_y=fix(n/2);  %获取圆心坐标
img=zeros(m,n);  %提前定义滤波后的频谱，提高运行速度
for j=1:n
    for i=1:m
        d=sqrt((i-O_x)^2+(j-O_y)^2);    %计算两点之间的距离
        H(i,j)=1/((1+r/d)^4);   %布特沃斯高通滤波器，这里设n=2，即2阶滤波器，n可调
        img(i,j)=H(i,j)*img_f(i,j); %滤波
    end
end
img=ifftshift(img);            %傅里叶反变换
img=uint8(real(ifft2(img)));   %取实数部分
figure
imshow(img)

---

3.高斯高通滤波

高斯高通滤波函数：
$$H(u,v)=1-e^{-D^2(u,v)/2D_0^2}$$

---

代码如下：

im1=imread('椒盐噪声图像.jpg');
im1=rgb2gray(im1);
figure
imshow(im1)

r=80;   %设置滤波圆半径参数为80，可调
img_f=fftshift(fft2(double(im1)));  %傅里叶变换得到频谱
[m,n]=size(img_f);
O_x=fix(m/2);  
O_y=fix(n/2);  %获取圆心坐标
img=zeros(m,n);  %提前定义滤波后的频谱，提高运行速度
for j=1:n
    for i=1:m
        d=sqrt((i-O_x)^2+(j-O_y)^2);    %计算两点之间的距离
         H(i,j)=1-exp((-d.^2)/(2*(r).^2));   %高斯高通滤波器
        img(i,j)=H(i,j)*img_f(i,j); %滤波
    end
end
img=ifftshift(img);            %傅里叶反变换
img=uint8(real(ifft2(img)));   %取实数部分
figure
imshow(img)

---

4.小结

频率域除了上述三种高通滤波外，还有拉普拉斯锐化增强、同态滤波等方法，将在后面给出。

---