MATLAB 实现图像的傅立叶变换

作者：timerring

2023-03-03
山东
本文字数：1821 字
阅读完需：约 6 分钟

文章和代码以及样例图片等相关资源，已经归档至【Github 仓库：digital-image-processing-matlab】或者公众号【AIShareLab】回复 数字图像处理 也可获取。

目的

1.掌握二维 DFT 变换及其物理意义

2.掌握二维 DFT 变换的 MATLAB 程序

3.空域滤波与频域滤波

原理

1.应用傅立叶变换进行图像处理

傅里叶变换是线性系统分析的一个有力工具，它能够定量地分析诸如数字化系统、采样点、电子放大器、卷积滤波器、噪音和显示点等的作用。通过实验培养这项技能，将有助于解决大多数图像处理问题。对任何想在工作中有效应用数字图像处理技术的人来说，把时间用在学习和掌握博里叶变换上是很有必要的。

2.傅立叶（Fourier）变换的定义

对于二维信号，二维 Fourier 变换定义为：

F (u, v) = \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} f (x, y) e^{- j 2 π (u x + v y)} d x d y

逆变换：

f (x, y) = \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} F (u d, v) e^{j 2 π (u x + v y)} d u d v

二维离散傅立叶变换为：

F (u, v) = x = 0 \sum M - 1 y = 0 \sum N - 1 f (x, y) e^{- j 2 π (\frac{u x}{M} + \frac{v y}{N})}

逆变换：

f (x, y) = \frac{1}{M N} u = 0 \sum M - 1 v = 0 \sum N - 1 F (u, v) e^{j 2 π (\frac{u x}{M} + \frac{v y}{N})}

图像的傅立叶变换与一维信号的傅立叶变换变换一样，有快速算法，具体参见参考书目，有关傅立叶变换的快速算法的程序不难找到。实际上，现在有实现傅立叶变换的芯片，可以实时实现傅立叶变换。

利用 MATLAB 实现数字图像的傅立叶变换

A．

I=imread('D:\pic\DIP3E_CH03\Fig0316(3)(third_from_top).tif'); %读入原图像文件imshow(I); %显示原图像fftI=fft2(I); %二维离散傅立叶变换sfftI=fftshift(fftI); %直流分量移到频谱中心RR=real(sfftI); %取傅立叶变换的实部II=imag(sfftI); %取傅立叶变换的虚部A=sqrt(RR.^2+II.^2);%计算频谱幅值A=(A-min(min(A)))/(max(max(A))-min(min(A)))*225; %归一化figure %设定窗口imshow(A); %显示原图像的频谱

复制代码

B．

傅立叶变换在图像处理，特别是在图像增强、复原和压缩中，扮演着非常重要的作用。实际中一般采用一种叫做快速傅立叶变换（FFT）的方法，MATLAB 中的 fft2 指令用于得到二维 FFT 的结果，ifft2 指令用于得到二维 FFT 逆变换的结果。

近似冲击函数的二维快速傅立叶变换（FFT）

x=1:99;y=1:99;[X,Y]=meshgrid(x,y);A=zeros(99,99);A(49:51,49:51)=1;B=fft2(A);subplot(1,2,1),imshow(A),xlabel('空域图象');subplot(1,2,2),imshow(B),xlabel('时域图象');figuresubplot(1,2,1),mesh(X,Y,A),xlabel('空域'),grid on;subplot(1,2,2),mesh(X,Y,abs(B)),xlabel('时域'),grid on;

复制代码

空域滤波与频域滤波

% 用于频域滤波的m函数function g=dftfilt(f,H)F=fft2(f,size(H,1),size(H,2));g=real(ifft2(H.*F));g=g(1:size(f,1),1:size(f,2));
%function PQ=paddedsize(AB,CD,PARAM)if nargin==1PQ=2*AB;elseif nargin ==2 & ~ischar(CD)PQ=AB+CD-1;PQ=2*ceil(PQ/2);elseif nargin == 2m=max(AB);P=2^nextpow2(2*m);PQ=[P,P];elseif nargin == 3m=max([AB CD]);P=2^nextpow2(2*m);PQ=[P,P];elseerror('Wrong number of inputs')end

%图像f的傅里叶频谱f=imread(' D:\pic\DIP3E_CH04\Fig0438(a)(bld_600by600).tif');F=fft2(f);S=fftshift(log(1+abs(F)));%S=gscale(S);% gscale 函数参照  function g=gscale(f,varargin)imshow(S)
%使用函数fspecial生成空间滤波器h=fspecial(‘sobel’)freqz2(h) %查看相应频域滤波器图形PQ=paddedsize(size(f));H=freqz2(h,PQ(1),PQ(2));H1=ifftshift(H);
imshow(abs(H),[])figure,imshow(abs(H1),[])
gs=imfilter(double(f),h);% 生成滤波后的图像，并默认采用了0进行边界填充gf=dftfilt(f,H1);
imshow(gs,[])figure,imshow(gf,[])
figure,imshow(abs(gs),[])figure,imshow(abs(gf),[])
%创建一幅阀值2值图像figure,imshow(abs(gs)>0.2*abs(max(gs(:))))figure,imshow(abs(gf)>0.2*abs(max(gf(:))))d=abs(gs-gf);max(d(:))min(d(:))