Matlab图像处理系列4———图像傅立叶变换与反变换
来源:互联网 发布:天下3 mac 编辑:程序博客网 时间:2024/06/05 01:14
注:本系列来自于图像处理课程实验,用Matlab实现最基本的图像处理算法
1.Fourier变换
(1)频域增强
除了在空间域内可以加工处理图像以外,我们还可以将图像变换到其他空间后进行处理,这些方法称为变换域方法,最常见的变换域是频域。
使用Fourier变换把图像从空间域变换到频域,在频域内做相应增强处理,再从频域变换到空间域得到处理后的图像。
我们这里主要学习Fourier变换和FFT变换的算法,没有学过通信原理,我对信号、时域分析也不是很清楚。
2.FFT算法
(1)离散Fourier变换,DFT
函数
很显然求出所有的长度为N的信号的DFT需要
(2)快速Fourier变换,FFT
快速傅立叶变换FFT是利用单位复数根的特殊性质(消去引理和折半引理,见《算法导论》(第3版中文版)P532详细论述),在
FFT有两种基本实现方式:
- 递归FFT
- 迭代FFT,也叫FFT蝶式运算
递归FFT由于递归栈开销大且容量有限等缺陷(但理解容易),一般计算采取迭代FFT实现。
(3)迭代FFT
直接给出算法导论版本的迭代FFT算法:
其中求逆序数拷贝的函数为:
FFT采用折半迭代的思想,因此速度能降低到
(4)迭代FFTMatlab实现
Matlab有fft函数,我们也可以自己实现:
function [ fft_m ] = IterativeFFT( vec ) clear i; n = length(vec); fft_m = BitReverseCopy(vec); for s = 1 : log2(n) m = power(2, s); wm = exp(- 2 * pi * i / m); for k = 0 : m : n - 1 w = 1; for j = 0 : m / 2 - 1 t = w * fft_m(k + j + m / 2 + 1); u = fft_m(k + j + 1); fft_m(k + j + 1) = u + t; fft_m(k + j + m / 2 + 1) = u - t; w = w * wm; end end endend
BitReverseCopy函数如下:
function [ copy ] = BitReverseCopy( vec ) n = length(vec); copy = zeros(1, n); bitn = log2(n); for i = 0 : n - 1 revi = bin2dec(fliplr(dec2bin(i, bitn))); copy(revi + 1) = vec(i + 1); endend
需要特别注意的是:
- 一般给出的FFT算法伪代码都是基于下标从零开始的数组,写在Matlab需要考虑映射关系
clear i
是为了怕之前有变量i和复数记号i混淆,清楚Matlab workspace中的缓存- 默认vec是double类型的!因为中间采用许多double类型运算
3.图像的二维Fourier变换
(1)二维DFT
二维DFT定义公式为:
计算一个频域点需要
(2)将二维DFT分解为两个一维DFT
Fourier变换的变换核(公式中和
先对二维矩阵的每一列做一维DFT:
再对变换后的矩阵的每一行做一维DFT:
最后以
(3)用一维FFT实现二维FFT
同样的我们可以用两个一维FFT实现二维FFT,时间复杂度
function [ mfft2 ] = JCGuoFFT2( data ) h = size(data, 1); w = size(data, 2); mfft2 = data; if power(2, log2(h)) ~= h || power(2, log2(w)) ~= w disp('JCGuoFFT2 exit: h and w must be the power of 2!') else for i = 1 : h mfft2(i, :) = IterativeFFT(mfft2(i, :)); end for j = 1 : w mfft2(:, j) = IterativeFFT(mfft2(:, j)); end endend
代码很简单,先做FFT行变换再做FFT列变换。之前忘记提到,我这里实现的都是长度必须是2的次方的Fourier变换,因此有时候会做一些长度规范检查。
(4)变换结果
经过JCGuoFFT2的二维傅立叶变换,我们可以得到复平面内各个点的复数值,那么显示在图像中需要先求出幅值:
pic1_fft = JCGuoFFT2(double(pic1));pic1_fft_amp = abs(pic1_fft);
在对幅值做一次log变换得到较好的频域图像:
pic1_fft_amp_log = log(1 + pic1_fft_amp);
绘制结果如下图:
(5)低频信号移到图像中心点
由于复数运算的周期特性,图像的Fourier变换在复平面上是完全对称的,可以想象平面是由无限多个上图(右)频域图像拼接而成的二维阵列。一般研究频域图像是把低频部分,也就是变换后的边缘部分移到图像中心点,Matlab提供fftshift
函数完成平移。
平移的思路有两个
- 通过Fourier变换平移定理先把原始图像做变换再做FFT
- 先做FFT后再根据频域图像的对称性做对称变换
查阅Matlab文档发现它是采用第二种方法,对图像做以下子矩阵交换:
那么我们可以很容易的写出自己的fftshift
:
function [ after ] = FFTShift( before ) h = size(before, 1); w = size(before, 2); after = before; if power(2, log2(h)) ~= h || power(2, log2(w)) ~= w disp('FFTShift exit: h and w must be the power of 2!') else hh = h / 2; hw = w / 2; after(1 : hh, 1 : hw) = before(hh + 1 : h, hw + 1 : w); after(1 : hh, hw + 1 : w) = before(hh + 1 : h, 1 : hw); after(hh + 1 : h, 1 : hw) = before(1 : hh, hw + 1 : w); after(hh + 1 : h, hw + 1 : w) = before(1 : hh, 1 : hw); endend
将低频部分平移到中心点后结果为:
4.图像的二维反Fourier变换
(1)一维逆DFT和一维逆FFT
一维离散傅立叶变换的逆变换是将e的指数部分变号,然后整体除以长度N(Fourier变换与逆变换关于符号、系数有很多种组合的定义,但他们都是等价且对称的。这里的定义配合Matlab做fft实际效果。):
同样我们可以根据iDFT的定义推演iFFT的算法,其迭代版本的Matlab实现如下:
function [ ifft_m ] = IterativeIFFT( vec ) clear i; n = length(vec); ifft_m = BitReverseCopy(vec); for s = 1 : log2(n) m = power(2, s); wm = exp(2 * pi * i / m); for k = 0 : m : n - 1 w = 1; for j = 0 : m / 2 - 1 t = w * ifft_m(k + j + m / 2 + 1); u = ifft_m(k + j + 1); ifft_m(k + j + 1) = u + t; ifft_m(k + j + m / 2 + 1) = u - t; w = w * wm; end end end ifft_m = ifft_m ./ n;end
(2)二维逆FFT
二维逆FFT和二维FFT的思路一致,对所有行和列分别作一次iFFT即可:
function [ mifft2 ] = JCGuoIFFT2( data ) h = size(data, 1); w = size(data, 2); mifft2 = data; if power(2, log2(h)) ~= h || power(2, log2(w)) ~= w disp('JCGuoIFFT2 exit: h and w must be the power of 2!') else for i = 1 : h mifft2(i, :) = IterativeIFFT(mifft2(i, :)); end for j = 1 : w mifft2(:, j) = IterativeIFFT(mifft2(:, j)); end endend
(3)逆FFT结果
对之前Rect1.bmp用JCGuoFFT2变换的到的Fourier变换(非shift和log之后、非仅幅度部分)作FFT2逆变换可以直接得到原图像:
这幅图有多个结果,可以看title知道每个结果的含义,图2-1是用JCGuoIFFT2做傅立叶反变换的结果,得到的图像和原图像、和Matlab ifft2反变换后的图像都是一致的。
5.幅频特性与相频特性
(1)对振幅和相位单独做逆FFT
如果我们只把图像Fourier变换的振幅部分和相位部分单独做二维逆FFT,可以直观的感受人眼对图像幅频特性和相频特性的敏感度。
复数
相位角和相位定义为:
对相位反变换需要在乘以一个系数(我是调出来的,针对图像,肯定可以自动的做均衡化):
pic2_fft_angle = angle(pic2_fft);clear i;tmp = 10000 * exp(i * pic2_fft_angle);pic2_ifft_angle = uint8(JCGuoIFFT2(tmp));
对振幅和相位单独做逆FFT结果如下:
(2)人眼敏感度
观察上图,幅频特性主要涵盖了图像颜色的分布,相频特性主要刻画了图像的边界信息。人眼对图像的相频特性更加敏感,看相频特性能够大概知道图像内容。
6.Fourier变换的旋转定理
(1)Fourier变换旋转定理
(2)结果
Rect2.bmp是Rect1.bmp旋转45度的示意图(注:原Rect2.bmp是二值的,做了预处理,但其本身边界不平滑,导致效果不太好,但不影响观察旋转定理):
我们可以看到幅度FFT正变换和相位FFT你变换都是旋转了45度,但是幅度FFT逆变换区别较大。
- Matlab图像处理系列4———图像傅立叶变换与反变换
- 傅立叶变换与MATLAB图像处理
- 简单图像处理——傅立叶变换
- 傅立叶变换与图像处理
- QT 实现图像处理-傅立叶变换、傅立叶反变换、平滑、锐化与模板匹配
- QT实现图像处理-傅立叶变换、傅立叶反变换、平滑、锐化与模板匹配
- 图像处理傅立叶变换
- C#|图像快速傅立叶变换与反变换
- 傅立叶变换与图像处理的关系
- 图像处理中的傅立叶变换
- 图像处理中的傅立叶变换
- 图像处理中的傅立叶变换
- 图像处理中的傅立叶变换
- 图像处理复习2——图像傅立叶变换和频域滤波
- 图像的灰度变换——图像旋转、图像的反色处理、对比度拉伸
- 图像傅立叶变换 概念
- 图像傅立叶变换
- 图像的傅立叶变换
- 《Unix编程艺术》读书笔记(1)
- DataStartSignal API
- 使用zendstudio10.6轻松创建符合官方推荐目录结构的zendframework2应用程序
- Python中步长索引解析
- CodeForces 550CDivisibility by Eight(暴力)
- Matlab图像处理系列4———图像傅立叶变换与反变换
- 工厂设计模式
- DataUtil API
- 基于XMPP协议的低传输负载的即时通信方法及其系统 -专利
- 让人头疼的XML文档
- 营销型网站对网站所提供的优化方案包括哪些?
- VS error C1083: 无法打开包括文件:“gl\glew.h”: No such file or directory\
- UVa 1632 Alibaba
- Android常用工具类