Matlab图像处理系列——频率域图像增强之快速傅里叶变换及实现

最新推荐文章于 2024-09-05 11:43:21 发布

小灰灰的FPGA

最新推荐文章于 2024-09-05 11:43:21 发布

阅读量2.6k

点赞数 27

分类专栏： Matlab图像仿真系列文章标签： matlab 图像处理 fpga开发

本文链接：https://blog.csdn.net/m0_50111463/article/details/134736870

版权

Matlab图像仿真系列专栏收录该内容

14 篇文章 14 订阅

订阅专栏

微信公众号上线，搜索公众号小灰灰的FPGA,关注可获取相关源码，定期更新有关FPGA的项目以及开源项目源码，包括但不限于各类检测芯片驱动、低速接口驱动、高速接口驱动、数据信号处理、图像处理以及AXI总线等
在这里插入图片描述
本节目录

一、快速傅里叶变换FFT
1、FFT变换的意义
2、FFT算法
3、按时间抽取的基-2FFT算法
二、快速傅里叶变换FFT的Matlab仿真
三、Matlab仿真源码

本节内容
一、快速傅里叶变换FFT
1、FFT变换的意义
快速傅里叶变换，是一种能够快速计算离散傅里叶变换的高效算法。
计算离散域上的傅里叶变换，DFT及反变换对：
在这里插入图片描述

计算上式中每个u值对应的F(u),需要N次复数乘法和N-1次复数加法，为了计算长度为N的序列的快速傅里叶变换，共需要执行N^2次复数乘法和N(N-1)次复数加法。实现一次复数相加至少需要执行两次实数相加，执行一次复数相乘可能需要至多4次实数乘法和2次实数加法。直接处理图像数据，运算量会大得惊人，无法实时处理。
变量WN具有周期性和对称性：
在这里插入图片描述

快速傅里叶变换的基本思路：通过将较长的序列转换成相对短得多的序列来减少运算量。利用W矩阵元素的周期性和对称性，将DFT运算中的某些项合并计算，计算半个W序列。将一个长度为N的序列分解成两个长度为N/2的序列并分别计算DFT。
2、FFT算法
FFT算法的基本思路分为两类：一是按照时间抽取的快速傅里叶算法，另一类是按照频率抽取的快速傅里叶算法。
按照时间抽取的快速傅里叶算法，基于将输入序列f(x)分解，抽取成较短的序列，然后从这些序列的DFT中求得输入序列的F(u)的方法。抽取后的较短序列仍然可以分解，最终只需要计算一个很短序列的DFT。需要注意序列的长度是2的整数次幂，高效地进行抽取和运算。
按照频率抽取的快速傅里叶算法，基于将输出序列F(u)分解，抽取成较短的序列，并从f(x)计算这些分解后序列的DFT。抽取后的较短序列仍然可以分解，最终只需要计算一个很短序列的DFT。需要注意序列的长度是2的整数次幂，高效地进行抽取和运算。
基-2FFT,序列长度是2的整数次幂时的DFT运算。
3、按时间抽取的基-2FFT算法
基-2FFT算法，N/2一般应是偶数，得到的序列可以继续分解，分解过程一直持续到每个序列只需要2点的DFT,只需计算这一DFT值，这一运算称为FFT的基本运算，即蝶形运算，蝶形运算的基础单元示意图：
在这里插入图片描述

对于基-2FFT，序列长度为N=2^L,且N为偶数。将序列按照项数的奇偶分成两组：
在这里插入图片描述

上式是一个递推公式，FFT蝶形运算的依据，一个偶数长度序列的傅里叶变换可以通过它的奇数项和偶数项的傅里叶变换得到，从而可以将输入序列分成两部分分别计算并按照公司相加/相减。
一个8点按时间抽取的FFT的算法：
在这里插入图片描述

一对复数进行完蝶形运算后，不需要再次保留输入的复数对，也就是输出对可以和输入对放在相同的存储单元中，只需要和输入序列大小相等的存储单元即可，也就是一种“原位”运算。
如果要使用这种“原位”运算，输入序列就必须按照倒序存储。由于f(x)是逐次抽取的，必须对原输入码列倒转位序，得到的次序相当于原序列编号的二进制码位导致。也就是，将原序列编号按照二进制表示，并且将二进制的所有位次序颠倒，就得到了在实际中输入序列中应该使用的排序位置。

在这里插入图片描述

二、快速傅里叶变换FFT的Matlab仿真
Matlab提供了fft2和ifft2函数，分别计算二维傅里叶变换和反变换。同时也提供了fftshift，将傅里叶频谱图中的零频点移动到频谱图的中心位置。
abs()可得到幅度谱，angle()可得到相位谱。
fft2函数输出的频谱分析数据，是按照原始计算所得到的顺序来排列频谱的，并没有以零频为中心排列，fftshift函数，利用频谱的周期性特点，将输出图像的一半平移到另一端，从而使零频被移动到图像的中间。
在这里插入图片描述

三、Matlab仿真源码

clear all;
clc;
close all;
%读取图像
P1=imread('people1.png');
P1=rgb2gray(P1);
P2=imread('people2.png');
P2=rgb2gray(P2);
%fft变换
P1f=fft2(P1);
P2f=fft2(P2)
%将零点移动到中心
%对幅度做对数变换以压缩动态范围
temp1=log(1+abs(fftshift(P1f)));
temp2=log(1+abs(fftshift(P2f)));
%傅里叶反变换
P11=abs(ifft2(P1f));
P22=abs(ifft2(P2f));
%交换相位谱
P12=abs(P1f).*cos(angle(P2f))+abs(P1f).*sin(angle(P2f)).*i;
P21=abs(P2f).*cos(angle(P1f))+abs(P2f).*sin(angle(P1f)).*i;
%傅里叶反变换
P12=abs(ifft2(P12));
P21=abs(ifft2(P21));
figure;
subplot(3,2,1);
imshow(P1);
title('people1原图');
subplot(3,2,2);
imshow(P2);
title('people2原图');
subplot(3,2,3);
imshow(temp1,[]);
title('FFT-P1');
subplot(3,2,4);
imshow(temp2,[]);
title('FFT-P2');
subplot(3,2,5);
imshow(P11,[]);
title('ifft-P11');
subplot(3,2,6);
imshow(P22,[]);
title('ifft-P22');
figure;
subplot(3,2,1);
imshow(P1);
title('people1原图');
subplot(3,2,2);
imshow(P2);
title('people2原图');
subplot(3,2,3);
imshow(temp1,[]);
title('FFT-P1');
subplot(3,2,4);
imshow(temp2,[]);
title('FFT-P2');
subplot(3,2,5);
imshow(P12,[]);
title('交换相位谱P12');
subplot(3,2,6);
imshow(P21,[]);
title('交换相位谱P21');