Matlab实现基于频域对二维信号的低通滤波

基于频域的低通滤波（二维信号——图像）

算法分析

傅里叶变换，将灰度图由f(x,y)->F(u,v)(空域转频域),得到图像在频域中的频谱（在频谱中低频信号分布在频谱的四个角落，其余部分为高频信号，这样的分布难以滤出高频或低频信号）
中心化，将频谱F(u,v)中心化，将低频点移到频谱中心（这样就可以通过设置一个截止频率D0，来过滤信号）
遍历频谱图，使用巴特沃兹低通滤波器和高斯高斯低通滤波器进行低通滤波，计算滤波器函数h(u,v)与F(u,v)的乘积G(u,v)，直到频谱图遍历完 巴特沃兹低通滤波器和高斯低通滤波器，计算公式
反中心化，对滤波后的频谱G(u,v)的低频点移回到频谱的四周
傅里叶反变换，将第4步的结果傅里叶反变换G(u,v)->g(x,y)，取g(x,y)实数部分作为最后滤波之后的结果，虚数部分是浮点运算存在的误差造成的（虚数部分的绝对值很小，可以忽略不计），得到滤波后的空域图像

伪代码

x = 读入灰度图，若是对彩色图处理，将彩色图先转为灰度图
y = 图像加噪
f = 将y的数据类型转换成double，便于运算
F = 二维傅里叶变换（这里可以直接调用fft2函数），得到图像在频域里的频谱
F1 = 对频谱你中心化，将低频频谱由四周转换到中心，距离中心越远的频率越高
[M, N] = 求得频谱的大小
n = 巴特沃斯低通滤波器的阶数
D0 = 设置截止频率
m = 频谱中心行坐标
n = 频谱中心列坐标
遍历频谱
        D = 计算频谱中当前点与频谱中心点的距离
        h1 = 根据步骤3的公式1计算巴特沃兹低通滤波器
        h2 = 根据步骤3的公式2计算高斯低通滤波器
        G1(i,j) = 计算巴特沃兹低通滤波器当前像素点灰度值的乘积
        G2(i,j) = 计算高斯低通滤波器当前像素点灰度值的乘积
G1 = 将巴特沃兹低通滤波后的频谱反中心化
G11 = 将频谱G1二维傅里叶反变换后，取实部
G2 = 将高斯低通滤波后的频谱反中心化
G22 = 将频谱G2二维傅里叶反变换后，取实部
%将巴特沃兹低通滤波后的空域图像G11和高斯低通滤波后的空域图像G22显示出来

代码

A = imread('lena.png');
x = rgb2gray(A);
y=imnoise(x,'gaussian',0,0.02);  %加入均值为0，方差为0.02的高斯噪声
f = double(y); %将数据类型转换成double
F = fft2(f); %二维傅里叶变换
F1 = fftshift(F); %中心化，将低频频谱由四周转换到中心，距离中心越远的频率越高
[M, N] = size(F);
n = 2;  %巴特沃斯低通滤波器的阶数
D0 = 50;  %设置截止频率为50
m = fix(M/2);
n = fix(N/2);
for i=1:M  
    for j=1:N
        D = sqrt((i-m)^2+(j-n)^2);
        h1 = 1 / (1 + (D / D0)^(2 * n));  %计算巴特沃兹低通滤波器
        h2 = exp(-(D.^2)./(2 * (D0^2)));  %计算高斯低通滤波器
        G1(i,j) = h1 * F1(i,j);  
        G2(i,j) = h2 * F1(i,j);
    end
end
G1 = ifftshift(G1);  %将巴特沃兹低通滤波后的频谱反中心化
G11 = uint8(real(ifft2(G1)));  %取二维傅里叶反变换后的实部
G2 = ifftshift(G2);  %将高斯低通滤波后的频谱反中心化
G22 = uint8(real(ifft2(G2)));  %取二维傅里叶反变换后的实部
subplot(2,4,1),imshow(x),title('原图');
subplot(2,4,2),imshow(y),title('加入高斯噪声后');
subplot(2,4,3),imshow(G11),title('二阶巴特沃斯低通滤波后图像');
subplot(2,4,4),imshow(G22),title('高斯低通滤波后图像');
subplot(2,4,5),imshow(log(1+abs(F)),[ ]),title('加噪图像的频谱');
subplot(2,4,6),imshow(log(1+abs(F1)),[ ]),title('频谱中心化');
subplot(2,4,7),imshow(log(1+abs(G1)),[ ]),title('巴特沃兹低通滤波器频谱反中心化');
subplot(2,4,8),imshow(log(1+abs(G2)),[ ]),title('高斯低通滤波器频谱反中心化');

实验结果

实验分析

由图所示，直接调用fft2函数将空域转到频域时，低频点分布在频谱的四周，不管是低通还是高通滤波都比较难限制，所以我们可以通过调用fftshift将频谱中心化，把低频信号移动到频谱中心，这样我们可以通过指定D0（截止频率），通过设置低通滤波器计算频谱上的一点到频谱中心的距离D(U,V) 与D0比较，大于D0舍弃，只保留小于D0的低频点
经过以上在频域中低通滤波，将频域反变换回空域时要将低频信号和高频信号反中心化回去，再进行反变换
原理总结：傅里叶变换后可以得到信号中有哪些频率的信号，将需要滤除的频率成分的幅值置为0，然后再经过傅里叶反变化就得到滤波的目的
为什么要将时域转换到频域？在时域中无法转换的数学函数，转换到频域能简单的实现