频率滤波（2）

function [f_filt]=freq_filter(f,h)f=mat2gray(f); % 归一化 转化为double h=double(h);pad_h=zeros(size(f)+size(h)-1);  % 在频率域滤波 空域的图像和滤波器都要填充，如果不填充会产生缠绕误差，填充后的大小pad_f=pad_h; % 填充后滤波器的大小与填充后图像的大小是相同的%但是填充零的话，相当于在空域乘以了一个距函数 即f(x,y)*rect(x,y)，那么在频域等价于%conv2(F(u,v),F(rect(x,y))[卷积定理]  众所周知 F（rect(x,y))  是一个sinc(x,y) 函数%sinc(x,y) 函数是向两边无限延伸的函数，造成有sinc 函数的高频分量导致的频率泄漏，泄漏会在图像上产生块效应%填充零是一种急剧过渡，后来人们就使用各种窗，来平滑过渡n=size(pad_h,1); %m=size(pad_h,2);nf1=fix((n-size(f,1))/2);mf1=fix((m-size(f,2))/2);nf2=fix((n+size(f,1))/2);mf2=fix((m+size(f,2))/2);np1=fix((n-size(h,1))/2);mp1=fix((m-size(h,2))/2);np2=fix((n+size(h,1))/2);mp2=fix((m+size(h,2))/2);pad_h(np1:np2-1,mp1:mp2-1)=pad_h(np1:np2-1,mp1:mp2-1)+h; % 填充算符（滤波器）pad_f(nf1:nf2-1,mf1:mf2-1)=pad_f(nf1:nf2-1,mf1:mf2-1)+f;% 图像填充f_filt=ifftshift(ifft2(fft2(pad_f).*fft2(pad_h)));  % 重建图像，并反变换到中心f_filt=real(f_filt); % 取实部% figure,imshow(real(f_filt),[]);end

clc;clear;f=imread('C:\Program Files\MATLAB\R2013a\bin\Original_Images\DIP3E_Original_Images_CH04\Fig0438(a)(bld_600by600).tif');% sobel=double([-1,0,1;-2,0,2;-1,0,1])';w1=double([1,2,1;0,0,0;-1,-2,-1]); %Sobel operator 3x3 convolution kernels% 水平 Sobel_H   kron([1,0,-1]',[1,2,1])w2=double([1,0,-1;2,0,-2;1,0,-1]); %Sobel operator 3x3 convolution kernels% 垂直 sobel_V   kron([1,2,1]',[1,0,-1])% sobel=double(fspecial('sobel'));  %水平 Sobel_H   kron([1,0,-1]',[1,2,1])f_filt1=freq_filter(f,w1);% sobel H operator f_filt2=freq_filter(f,w2);% vfigure,imshow(f,[]);title('ori image')figuresubplot(2,2,1),imshow(f_filt1,[]);title('sobel H 算子频域滤波')subplot(2,2,2),imshow(f_filt2,[]);title('sobel V 算子频域滤波')f_filt3=conv2(f,w1);   % sobel H算子空域滤波f_filt4=conv2(f,w2);   % sobel V算子空域滤波subplot(2,2,3),imshow(real(f_filt3),[]);title('sobel H 算子空域滤波')subplot(2,2,4),imshow(real(f_filt4),[]);title('sobel V 算子空域滤波')