MATLAB频域处理-傅里叶变换和滤波

主要demo：

二维傅立叶变换，二维FFT，频域低通滤波，高通滤波，拉普拉斯算子

---

      f=imread('你的图片');
      F=fft2(f);            %对图像f进行傅里叶变换
      S=abs(F);            %S是F的频谱
      imshow(S,[]);       %显示频谱，能量集中子在四个角
      FC=fftshift(F);       %频谱移动到中心
      figure,imshow(FC,[]);        %显示频谱，能量移动到中心
      S2=log(1+abs(FC));    %对频谱作对数变换，拓展其动态范围
      figure,imshow(S2,[]);        %显示变换后的频谱，频谱的细节部分也能清楚显示

f =imread('你的图'); 
imshow(f);
[M,N]=size(f);   % M,N分别是图像f的宽度和高度

%对f进行傅立叶变换，没有对f进行扩充填0操作
F=fft2(f);
sig=10;         % sig是高斯滤波器的方差

% lpfilter 产生一个频域低通滤波器
H =lpfilter('gaussian', M, N, sig);
G=H.*F;    % 频域滤波，将滤波器和F相乘
% 滤波结果进行傅立叶逆变换并取其实部，得到空间域的滤波结果
g =real(ifft2(G));
% 由于折叠效应，g的图像水平边缘有些模糊，但是垂直边缘并不模糊
figure, imshow (g,[]);
PQ=paddedsize (size(f));  %PQ是扩充后的图像尺寸，默认为（M*2, N*2）
Fp=fft2(f, PQ(1), PQ(2));   %傅立叶变换并扩充到指定尺寸
Hp=lpfilter('gaussian', PQ(1),PQ(2), 2*sig);  %Hp是指定尺寸的低通滤波器
Gp=Hp.*Fp;     %频域低通滤波
gp=real(ifft2(Gp));  % 傅立叶逆变换，得到滤波结果
figure, imshow(gp,[  ]);  % 傅立叶逆变换后的结果，原图像只占据左上角1/4
gpc=gp(1:size(f,1), 1:size(f,2));  %只取逆变换的左上角部分
figure, imshow(gpc, [  ]);  %可以看到图像的垂直和水平边缘都变模糊

f =imread('你的图');
imshow(f);
F=fft2(f);       %对f做傅立叶变频
% S是对数变换并移到中心的傅立频谱
S=fftshift(log(1+abs(F)));
S=gscale(S);      %将S归一化到[0，1]之间
figure, imshow(S,[ ]);  %显示傅立频谱

h=fspecial ('sobel')  %是3*3 垂直方向Sobel滤波器
freqz2(h)     %显示与h相应的频域滤波器的图形
PQ= paddedsize(size(f));  %填充大小
H=freqz2(h, PQ(1), PQ(2)); %生成与h相应的频域滤波器，并扩充到指定大小
H1=ifftshift(H);           %解除H的居中
figure, imshow(abs(H),[  ]);  %显示频域滤波器的频谱，能量集中在中心
figure, imshow(abs(H1),[  ]);  %解除居中操作后的频谱，能量集中在四个角上

gs=imfilter(double(f), h);    % 对f进行空间滤波操作，采取0填充边界
gf=dftfilt(f, H1);           %对f 进行频域滤波，dftfilt 函数见第4.6节
figure, imshow(gs, [ ]);   %显示空间滤波结果
figure, imshow(gf, [ ]);    %显示频域滤波结果

figure, imshow(abs(gs), [ ]);  %取gs的绝对值，消除负值
figure, imshow(abs(gs), [ ]);  %取gf的绝对值，消除负值

d=abs(gs-gf);     %对空间滤波结果和频域滤波结果取差值
max(d(:))               %打印差值的最大值
min(d(:))              %打印差值的最小值

f =imread('你的图');
imshow(f);
PQ=paddedsize(size(f));          %计算填充大小
%利用函数 dftuv 得到频域坐标网格
[U, V]=dftuv(PQ(1), PQ(2));
DO=0.05*PQ(2);               % DO是截至频率
F = fft2(f,PQ(1), PQ(2));        % 傅立叶变换
H=exp(-(U.^2+V.^2)/(2*(DO^2)));   % H是高斯低通滤波器
g =dftfilt(f,H);              % 高斯低通滤波的结果
figure, imshow(g,[]);
figure, imshow(fftshift (H), [ ] );  % 查看居中之后的低通滤波器
figure, imshow(log(1+abs(fftshift (F))), [ ] );  % 显示对数变换后的频谱
figure,mesh(fftshift(H(1:10:PQ(1), 1:10:PQ(2))));      % 以网格的形式显示低通滤波器
figure,surf(fftshift(H(1:10:PQ(1), 1:10:PQ(2))));       % 以表面图的形式显示低通滤波器

f =imread('你的图');
H=fftshift(hpfilter('ideal', 500, 500, 50));  % H是居中的高通滤波器
figure,mesh(H(1:10:500, 1 :10:500));          % 以网格的形式显示高通滤波器
figure,imshow(H, [ ]);                    % 以图像形式显示滤波器

PQ = paddedsize(size(f));  %填充大小
D0 = 0.05*PQ(1);          %截止频率
H = hpfilter('gaussian', PQ(1),PQ(2),D0);%高通滤波器
g = dftfilt(f, H);         %高通滤波
figure,imshow(g,[  ]);
%**********高频强调滤波***********************
f = imread('你的图');
figure,imshow(f);
PQ = paddedsize(size(f));
DO = 0.05 * PQ(1);
% HBW是巴特沃斯高通滤波器
HBW = hpfilter('btw',  PQ(1),  PQ(2),  DO,2);
H = 0.5+2 * HBW;           %H是高频强调滤波器
gbw = dftfilt(f,HBW);     %高通滤波
gbw=gscale(gbw);         %将值归一化到[o,1]
figure,imshow(gbw,[  ])  %显示高通滤波图像
ghf = dftfilt(f,H);          %高频强调滤波
ghf = gscale(ghf);           %值归一化
figure,imshow(ghf,[  ])   %显示高频强调滤波后的图像
ghe = histeq(ghf,256);     %直方图均衡化
figure,imshow(ghe,[  ])  %直方图均衡化后的图像