一种基于共轭对称性的二维二维ff正反变换实时处理方法

一种基于共轭对称性的二维二维ff正反变换实时处理方法

二维快速傅里叶变换（fft）是数据处理领域中常见的分析和计算工具。它在理论研究和工程应用中具有重要的价值。图像低频分析、频率滤波、卷积计算和相关计算等通常是典型的。但高分辨率图像二维快速傅里叶变换的数值计算量大且动态范围宽、输入/输出数据依赖性强、硬件实现复杂度高等特点,制约了其在图像处理领域的实时性应用。通常,利用多维傅里叶变换的可分性,图像的二维FFT可通过分别计算行方向和列方向的一维FFT来实现。应用相对成熟的一维FFT优化方法,上述行列算法在很大程度上可加快二维FFT的计算速度,但仍难以满足目前高分辨率图像二维FFT的实时处理要求。在实信号FFT中引入傅里叶变换的奇偶对称性,可减少二维FFT计算和内存使用,极具工程价值。另外,对于二维FFT的高速实现,通常以高速FPGA、高性能DSP、通用GPU芯片为核心计算器件的三种常用实现平台,并在系统设计中采用乒乓缓存、流水线、多核计算等并行处理技术。本文总结了常规的二维FFT行列算法和实现中的技术难点,并基于实际工程中二维图像均为实值信号的物理事实和实信号傅里叶变换的周期对称性和数据共轭对称性,给出高分辨图像二维FFT正/反变换处理方法,目的在于大幅缩减FFT计算量和数据存储需求。1维fft的计算所面临的挑战由一般的行列算法,可以得到常规的二维FFT正/反变换处理流程,并用图1中的模块化形式给出。其中模块a、b完成二维FFT正变换功能,模块c、d则完成二维FFT反变换,以二维IFFT表示。二维FFT的可分性没有计算顺序要求,先行后列或者先列后行的处理顺序不会影响正变换或反变换的输出结果,故模块a和b、模块c和d可以互换,只需保证它们各完成一次行和列方向的一维FFT或IFFT即可。图1中正变换采用了先行后列,而反变换采用了先列后行的顺序,这样有利于减少硬件中数据的寻址时间。二维FFT正/反变换的完成取决于一维FFT的处理速度,而大点数的一维FFT本身就是一个计算密集型问题。例如对N×N像素的图像,完成一次二维FFT正变换需要N2log2N次复数乘法和2N2log2N次复数加法。另一方面,一维FFT的计算是对所有点的输入数据进行乘积累加完成的,其本身有着很强的数据依赖性。而将二维FFT分解为行和列方向的一维FFT后,更加剧了这种数据依赖性。这在很大程度上增加了处理时间和硬件设计难度,不仅限制了数据的流水线处理,还需要在行、列方向一维FFT之间增加大容量存储,用来缓存一帧的复值图像数据。此外,长序列一维FFT的输出数据具有较宽的动态范围,主要与输入数据位宽和变换点数有关。因此在保证数据最高有效位的前提下,若使用定点数制计算高分辨图像的二维FFT,则精度损失较大,难以满足误差限要求。使用浮点数制可改善固定位宽对计算精度的影响,但代价是,在整个二维FFT处理过程中需要使用更长处理时延的浮点数乘加运算和更大的存储空间,故在硬件设计中需要对系统的数据格式的选择进行权衡。为保证变换结果的数值精度,本文的理论分析和硬件设计均基于IEEE-754标准的32位单精度浮点数格式。如上所述,一维FFT的使用令二维FFT的实现成为可能,但输入数据的规模性增长,使得一维FFT处理中的技术难点也被放大性地引入到了二维FFT处理中。概括起来,包括以下三个方面:a)大量的一维FFT计算使得二维FFT成为典型的计算密集型问题;b)傅里叶变换对全域数据的依赖性,导致了二维FFT行列变换间的处理时延和大容量缓存;c)为保证计算精度,浮点数的使用进一步加剧了计算和存储需求。本文正是针对上述三个主要的技术难点,在处理方法和硬件设计两个层面对高分辨图像的二维FFT实时处理进行改进和优化。2改进的基本原则2.1复序列的同化作用计算一维实信号的FFT,可以利用傅里叶变换的周期对称性将N点的实序列奇偶分离,构造一个N/2点的复序列,把N点的实序列FFT计算转换为计算一个N/2点的复序列FFT;然后将该复序列的N/2点FFT输出进行一定的运算与组合,获得原N点实序列的FFT输出。使用该方法可以缩减一维实序列FFT近一半的计算量。基本原理如下:将一个N点的序列x(n)按奇偶序号分为奇数序列g(n)和偶数序列h(n),其离散傅里叶变换分别记为G(k)和H(k)。根据FFT时间抽取基2算法:X(k)=G(k)+WkNH(k)k=0⋯N2−1(1)X(k)=G(k)+WΝkΗ(k)k=0⋯Ν2-1(1)X(k+N2)=G(k)−WkNH(k)k=0⋯N2−1(2)X(k+Ν2)=G(k)-WΝkΗ(k)k=0⋯Ν2-1(2)构造y(n)=g(n)+j·h(n),则y(n)为N/2点的复序列。y(n)的离散傅里叶变换Y(k)为Y(k)=∑n=0N/2−1(g(n)+j⋅h(n))⋅WnkN/2k=0⋯N2−1(3)Y(k)=∑n=0Ν/2-1(g(n)+j⋅h(n))⋅WΝ/2nkk=0⋯Ν2-1(3)令k=N2−kk=Ν2-k,代入式(3)可得Y(N2−k)=∑n=0N/2−1(g(n)+j⋅h(n))⋅Wn(N2−k)N/2=∑n=0N/2−1(g(n)+j⋅h(n))⋅W−nkN/2(4)Y(Ν2-k)=∑n=0Ν/2-1(g(n)+j⋅h(n))⋅WΝ/2n(Ν2-k)=∑n=0Ν/2-1(g(n)+j⋅h(n))⋅WΝ/2-nk(4)再令ω=2πN/2nkω=2πΝ/2nk,代入式(3)(4)且两式等号左右两边相加,可得Y(k)+Y(N2−k)=2⋅∑n=0N/2−1g(n)⋅cosω+2j⋅∑n=0N/2−1h(n)⋅cosω(5)Y(k)+Y(Ν2-k)=2⋅∑n=0Ν/2-1g(n)⋅cosω+2j⋅∑n=0Ν/2-1h(n)⋅cosω(5)将复序列Y(k)、G(k)、H(k)的实虚部分开,分别记为Yr(k)、Gr(k)、Hr(k)和Yi(k)、Gi(k)、Hi(k),代入式(5),则有Gr(k)=(Yr(k)+Yr(N2−k))/2(6)Gr(k)=(Yr(k)+Yr(Ν2-k))/2(6)Hr(k)=(Yi(k)+Yi(N2−k))/2(7)Ηr(k)=(Yi(k)+Yi(Ν2-k))/2(7)类似方法可得Gi(k)=(Yi(k)−Yi(N2−k))/2(8)Gi(k)=(Yi(k)-Yi(Ν2-k))/2(8)Hi(k)=−(Yr(k)−Yr(N2−k))/2(9)Ηi(k)=-(Yr(k)-Yr(Ν2-k))/2(9)因此,只要求出N/2点复序列y(n)的傅里叶变换Y(k),利用上述四个关系式(6)～(9)分离出g(n)和h(n)的傅里叶变换G(k)、H(k);再代入式(1)和(2),即可得到N点实序列x(n)的傅里叶变换X(k)。图1中的正变换模块a,其输入数据为实信号,因此可以利用该方法缩减一维FFT计算中近一半的计算量。而对于高分辨率图像而言,即N值较大的情况下,其节约的计算量是非常可观的。2.2维fft变换结果的深度学习实信号的一维FFT和二维FFT变换结果均具有中心共轭对称性。以一维FFT为例,如果信号f(n)是一个实信号,则可以得到F(k)=DFT{f(n)}=∑n=0N−1f(n)cos(2πNnk)−j⋅∑n=0N−1f(n)sin(2πNnk){f(n)}=∑n=0Ν-1f(n)cos(2πΝnk)-j⋅∑n=0Ν-1f(n)sin(2πΝnk)由上式可知,F(k)的实部是k的偶函数,虚部是k的奇函数,或者F(k)对k呈共轭对称性。这种对称性是由F(k)实虚部中正余弦函数的奇偶性决定的,而实际上k的取值范围是从0～N-1,这时对称中心位置应为N2Ν2,可表示为F(k)=F*(N-k)0≤k≤N-1(10)同样,这种中心共轭对称性扩展到二维FFT变换结果F(u,v)中,则如式(11)所示,二维数据的对称中心为(M2,N2)(Μ2,Ν2)。F(u,v)=F*(M-u,N-v)0≤u,v≤N-1(11)根据中心共轭对称性可知,图1中二维FFT正变换计算的中间结果F(x,v)和最终结果F(u,v)中存在着数据冗余。对于模块a的输出F(x,v),只需计算并存储其中半帧数据,其余数据则可以通过求取复数共轭来恢复,图2表示出了中间结果F(x,v)的这种对称性。同样道理,模块b的输出F(u,v),也只有半帧的数据具有独立意义。图3表示出了最终结果F(u,v)的这种对称性。由于中间结果F(x,v)与最终结果F(u,v)的数据依赖性只存在于每列数据之间,基于前面分析,如果只存储模块a的输出F(x,v)的前N/2列数据,那么就可以用这部分数据来计算出F(u,v)前N/2列数据,然后根据中心共轭对称性恢复出F(u,v)中的剩余数据,从而得到完整的二维FFT变换结果。这样处理的好处在于,利用实信号傅里叶变换的共轭对称性去除了二维FFT变换中间结果和最终结果中的冗余数据,使得模块a和b的输出数据存储量、模块b中一维FFT的计算量均减少了近一半。2.3中心共保持对称实值输入图像二维FFT正变换结果F(u,v)具有中心共轭对称性,由模块e对其进行频域处理。而图像频域处理中的很多计算操作,可使处理后输出F′(u,v)继续保持共轭对称性。例如,多数频域卷积或滤波模板也具有中心对称性,使用这些模板对F(u,v)进行频域处理后得到的F′(u,v)仍具有中心共轭对称性;又如,计算空域卷积的两幅实值图像,其二维FFT正变换结果均关于中心共轭对称,而两者做点乘后的计算结果也具有中心共轭对称性。因此,若频域处理结果可继续保持中心共轭对称性,则在模块e中,只需计算F(u,v)中前半帧或后半帧数据即可得到F′(u,v)。由于F′(u,v)具有共轭对称性,输入模块c进行列方向的一维FFT变换后,输出数据F′(x,v)同样关于中心共轭对称,进一步利用该性质又可以节省模块c中一半的一维FFT计算量和存储需求。模块d中计算也可以只利用F′(x,v)中一半的独立数据来完成,其基本原理是利用一维FFT周期对称性对2.2节中的过程进行反推,将原本N点的复数FFT转换为计算N/2点的复数FFT,从而使模块d中一维FFT的计算减少近一半。3当前处理方法和硬件实现3.1维fft的计算和存储应用上述方法对常规处理方法进行改进,可以得到二维FFT正/反变换实时处理方法,如图4所示。其处理流程如下所述:a)数据预处理I。输入大小为M×N的实值图像f(x,y),各行数据按奇偶序号分离;将各行奇、偶序列组成每行N/2点的复值数据fc(x,y)。b)模块a。对fc(x,y)各行数据进行N/2点的一维FFT计算,再利用式(3)(4)和(8)～(11)计算出每行的一维FFT输出序列;保存每行输出序列中前一半的数据,得到中间结果Fc(x,v)。c)模块b。对Fc(x,v)中每一列复值数据计算一维FFT,得到二维FFT变换最终结果中的前半帧数据Fc(u,v)。d)数据恢复Ⅰ。若频域处理可保持数据共轭对称性,则该模块直接输出Fc(u,v);否则,对Fc(u,v)中相应点求取复共轭,利用共轭对称性恢复完整的频域数据F(u,v)并输出。e)模块e。对Fc(u,v)或Fc(x,y)进行频域处理,得到反变换输入数据F′c(u,v)或F′(u,v)。f)模块c。对F′c(u,v)或F′(u,v)的每一列数据计算一维IFFT,得到反变换中间结果F′c(x,v)或F′(x,v)。g)数据预处理Ⅱ。若频域处理可保持数据共轭对称性,则利用F′c(x,v)计算出每行N/2点的复值序列Y(k)组成数据块F′d(x,v);否则,直接输出F′(x,v)。h)模块d。对输入数据F′d(x,v)或F′(x,v)中每一行计算一维IFFT,得到反变换最终结果f′d(x,v)或f′c(x,v),输出数据为复值。i)数据恢复Ⅱ。若频域处理可保持数据共轭对称性,则对f′d(x,v)各行数据进行虚、实部分离,组合成每行N点的实值图像f′(x,y);否则,对f′c(x,v)中各点取实部,得到实值图像f′(x,y)。相比常规处理方法,改进后方法只增加了四个功能模块,即数据预处理I和Ⅱ、数据恢复I和Ⅱ。而四个模块中的处理仅包含数据转存、改变数值符号和部分乘加运算,在硬件设计中均实现简单且耗时较少。只需引入较少的硬件成本,改进后的处理方法就可以缩减二维FFT正变换近一半的一维FFT计算和存储需求。在能够继续保持数据共轭对称性的情况下,又可以进一步缩减模块e、c、d中一半一维FFT计算量和存储。同时,为保证整个方法硬件实现的低复杂度和一维FFT的通用性,对模块a、b、c、d的功能无任何改动要求。3.2fpga实时实现和图像处理在dsp中的应用本文在改进的二维FFT正/反变换实时处理框架基础上,采用以支持浮点操作的DSP芯片ADSP-TS201为计算核心、共搭载4片TS201的图像处理平台,对高分辨率图像二维FFT正/反变换的实时实现和图像处理应用进行了设计开发。在不使用外部存储的情况下,使最大尺寸为512×512像素的图像处理得以在DSP内部实现,并且只以较少的数据存取时间作为代价(不到总处理时间的1%)。另外在系统设计中,将同一图像序列分为四路输入,使用4片TS201同时各自处理一路输入数据,再令四路处理数据同步输出恢复为一路图像输出序列,从而在实现图像序列的任务级并行处理,进一步提高系统的处理帧频。4图像性能对比为验证改进后实时处理方法的数据处理效率和硬件实现可行性,在上述DSP处理平台上进行硬件设计,并完成了以下两组实验。a)实现常规处理和实时处理两种方法下的二维FFT正/反变换。首先对不同分辨率的输入图像作二维FFT正变换,测算完成正变换的单片DSP耗时和四片DSP处理帧频,实验数据如表1所示;然后直接对正变换结果作二维FFT反变换,测算完成整个正/反变换过程的单片DSP耗时和四片DSP处理帧频,如表2所示。对比本次实验数据可知,对512×512像素的图像,单片DSP正变换处理时间仅24.3ms,一次正/反变换的处理