由卷积定理可知所有频域的滤波理论上都可以转化为空域的卷积操作。

给定频率域滤波器，可对其进行傅里叶逆变换得到对应的空域滤波器；滤波在频域更为直观，但空域适合使用更

小的滤波模板以提高滤波速度。因为相同尺寸下，频域滤波器效率高于空域滤波器，故空域滤波需要一个更小尺

寸的模板近似得到需要的滤波结果。

空域卷积

将模板在图像中逐像素移动，将卷积核的每个元素分别和图像矩阵对应位置元素相乘并将结果累加，累加和作为

模板中心对应像素点的卷积结果。通俗的讲，卷积就是对整幅图像进行加权平均的过程，每一个像素点的值，都

由其本身和邻域内的其他像素值经过加权平均后得到。在像素的处理上，是先将结果暂存在于一个副本，最后统

一拷贝，故不会出现处理顺序不同而结果不同的情况。

二维连续卷积的数学定义：

频率域是由傅里叶变换和频率变量 (u,v)定义的空间，频域滤波处理过程：先对图像进行傅里叶变换，转换至频率

域，在频域使用滤波函数进行滤波，最后将结果反变换至空间域。即：

高斯函数的特殊性：高斯函数傅里叶变换仍是高斯函数，但标准差已经变化，频域标准差越大（高斯函数越

宽），变换后空域标准差越小（高斯函数越窄）。

如频域函数对应的空域函数为

前面的系数在实际应用中常被忽略，因为是离散取样，不能使取样和为1，最后还要做归一化操作。

程序：

function filt=mygaussian(varargin)
%参数初始化，使用varargin处理可变参数情况
siz=varargin{1};%模板尺寸
if(numel(siz)==1)
    siz=[siz,siz];
end
std=varargin{2};%方差
centa = (siz(1)+1)/2;%此处不要取整
centb = (siz(1)+1)/2;

filt = zeros(siz(1),siz(2));
summ=0;
for i=1:siz(1)
    for j=1:siz(2)
        radius = ((i-centa)^2+(j-centb)^2);
        filt(i,j) = exp(-(radius/(2*std^2)));
        summ=summ+filt(i,j);
    end
end
filt=filt/summ;%归一化

测试：

执行mygaussian(4,1)得：

0.0181   0.0492    0.0492    0.0181

0.0492   0.1336    0.1336    0.0492

0.0492   0.1336    0.1336    0.0492

0.0181   0.0492    0.0492    0.0181

执行fspecial(‘gaussian‘,4,1)得：

0.0181   0.0492    0.0492    0.0181

0.0492   0.1336    0.1336    0.0492

0.0492   0.1336    0.1336    0.0492

0.0181   0.0492    0.0492    0.0181

查看fspecial的Matlab源码，将gaussian部分提取出来：

function h = hisfspecial(varargin)
siz=varargin{1};
if(numel(siz)==1)
    siz=[siz,siz];
end
std   = varargin{2};
siz   = (siz-1)/2;
[x,y] = meshgrid(-siz(2):siz(2),-siz(1):siz(1));
arg   = -(x.*x + y.*y)/(2*std*std);
h     = exp(arg);
h(h<eps*max(h(:))) = 0;
sumh = sum(h(:));
if sumh ~= 0,
    h  = h/sumh;
end;

与mygaussian类似，只是更具Matlab特色。

模板与图像卷积

直接利用Matlab提供的fspecial、imfilter。

I=imread('lena.bmp');
h=fspecial('gaussian',7,8);
O=imfilter(I,h,'replicate','same','conv');
subplot(1,2,1);imshow(I,[]);title('input');
subplot(1,2,2);imshow(O,[]);title('output'); 

具体的实现细节还有很多值得研究的，如边界处理、将二维高斯模板分解为两个一维模板等。可参考文后链接。

结果：

频域高斯低通滤波

在频率域，高斯函数表示为:

或写成

程序：

I=imread('lena.bmp');
I=double(I);
[M,N]=size(I);
m=(M+1)/2;
n=(N+1)/2; 
d0=30;                       %截止频率
F=fftshift(fft2(I));               %空域转频域，平移中心
H=zeros(M,N);
G=zeros(M,N);
for i=1:M
   for j=1:N
       d=(i-m)^2+(j-n)^2;
       H(i,j)=exp(-d/(2*d0^2));
       G(i,j)=H(i,j)*F(i,j);
    end
end
O=ifft2(ifftshift(G));                %平移中心，频域转空域
O=real(O);                        %取实数部分

subplot(1,2,1);imshow(I,[]);title('inpit');
subplot(1,2,2);imshow(O,[]);title('output');