《北京理工大学学报》：稀疏傅里叶变换理论及研究进展（仲顺安、王雄、王卫江、刘箭言）

第37卷第2期 北京理工大学学报 Vol 37 No. 2 017年2月 Transactions of Beijing Institute of Technology Feb.2017 稀疏傅里叶变换理论及研究进展 仲顺安1,王雄12,王卫江1,刘箭言1 (1.北京理工大学信息与电子学院,北京100081;2.重庆通信学院信息工程系,重庆400035) 摘要:稀疏傅里叶变换( sparse Fourier transform,SFT)是一种稀疏信号离散傅里叶变换的新算法,比传统快速 傅里叶变换( fast Fourier transform,FFT)更加高效.综述了SFT的理论框架、约束条件及频谱重排、窗函数滤波、 降采样FFT等关键技术问题,结合算法最新理论成果,归纳出4种不同的重构方法:哈希映射法、混叠同余法、相位 解码法、二分查找法,最后介绍了SFT理论的应用成果,并展望了其未来可能的发展方向 关键词:稀疏傅里叶变换;频谱重排;平坦窗函数;降采样FFT;哈希映射 中图分类号:TN911 文献标志码:A文章编号:1001-0645(2017)02-0111-08 DOI:10.15918/j, title001-0645.2017.02.001 Recent Advances in the sparse fourier Transform ZHONG Shun-an', WANG Xiong, 2, WANG Wei-jiang, LIU Jian-yan School of Information and Electronics, Beijing Institute of Technology, Beijing 100081, China; 2 Department of Information Engineering, Chongqing Communication Institute, Chongqing 400035, China) Abstract: Sparse Fourier transform (SFT) is a novel algorithm for discreting Fourier transform (DFT) on sparse signals, and is more efficient than the traditional fast Fourier transform (FFT). Reviewing the theoretical framework, restrictions and the key technical problems such as random spectrum permutation, window filtering and subsampled FFT, our different kinds of reconstruction means: hash mapping, aliasing-based search, phase decoding, binary search were introduced based on the latest theoretical achievements of the algorithms. Finally, some applications based on SFt were introduced, and its outlooks were presented. Key words: SFT; spectrum permutation; flat window function; subsampled FFT: Hash function 离散傅里叶变换( discrete Fourier transform,人们开始将目光集中在信号自身特征上来,研究发 DFT)是信号处理中最基础的理论算法之一口,计现,大量信号具有“稀疏性”:频域只有小部分大值 算DFT最广泛的方法是1965年由 Cooley等口提点,其余大部分点的值趋近于0.这一特征涵盖众多 出的快速傅里叶变换,其时间复杂度为O(NlbN).应用领域:图像和语音信号[、机器学习理论、布 与直接计算DFT相比,FFT大大简化了运算过程,尔函数分析、压缩感知-、宽带频谱感知 可将运算量缩短1~2个数量级,推动了信号处理技等.对于这种广泛存在且结构特殊的信号,能否实 术的革命性进展,被誉为二十世纪最具影响力的十现算法性能的突破呢? 大算法之一[3.半个世纪以来,研究人员不断追求 首先进行这一研究的是1993年文献[12]中设 FFT算法的改进并提出了若干变体,如原地FFT计的“哈达玛变换”快速算法,不久之后又出现了 算法和分裂基FFT算法[等,这些变体在某些硬件系列复数域上的研究[13-1.这些算法在一定的稀疏 实现上确有助益,但是复杂度并未得到明显改观.度范围内,性能超越FFT,其时间复杂度大都是由 收稿日期:2014-11-06 乍者简介:仲顺安(1957—),男,博士,教授,Emal: zhongsi@abit.edu.cm 通信作者:王卫江(1976—),男,博士,副教授, E-mail: wangweijiang@bit.edu.cr

收稿日期:2014 11 06 作者简介:仲顺安(1957—),男,博士,教授,E-mail:zhongsa@bit.edu.cn. 通信作者:王卫江(1976—),男,博士,副教授,E-mail:wangweijiang@bit.edu.cn. 第37卷 第2期 2017年2月 北 京 理 工 大 学 学 报 TransactionsofBeijingInstituteofTechnology Vol.37 No.2 Feb.2017 稀疏傅里叶变换理论及研究进展 仲顺安1, 王雄1,2, 王卫江1, 刘箭言1 (1. 北京理工大学 信息与电子学院,北京 100081;2. 重庆通信学院 信息工程系,重庆 400035) 摘 要:稀疏傅里叶变换(sparseFouriertransform,SFT)是一种稀疏信号离散傅里叶变换的新算法,比传统快速 傅里叶变换(fastFouriertransform,FFT)更加高效.综述了SFT的理论框架、约束条件及频谱重排、窗函数滤波、 降采样FFT等关键技术问题,结合算法最新理论成果,归纳出4种不同的重构方法:哈希映射法、混叠同余法、相位 解码法、二分查找法.最后介绍了SFT理论的应用成果,并展望了其未来可能的发展方向. 关键词:稀疏傅里叶变换;频谱重排;平坦窗函数;降采样 FFT;哈希映射 中图分类号:TN911 文献标志码:A 文章编号:1001-0645(2017)02-0111-08 DOI:10.15918/j.tbit1001-0645.2017.02.001 RecentAdvancesintheSparseFourierTransform ZHONGShun-an1, WANGXiong 1,2, WANG Wei-jiang 1, LIUJian-yan1 (1.SchoolofInformationandElectronics,BeijingInstituteofTechnology,Beijing100081,China; 2.DepartmentofInformationEngineering,ChongqingCommunicationInstitute,Chongqing400035,China) Abstract:SparseFouriertransform (SFT)isanovelalgorithmfordiscretingFouriertransform (DFT)onsparsesignals,andis moreefficientthanthetraditionalfastFouriertransform (FFT).Reviewingthetheoreticalframework,restrictionsandthekeytechnicalproblemssuchas randomspectrum permutation,windowfilteringandsubsampledFFT,ourdifferentkindsof reconstructionmeans:hashmapping,aliasing-basedsearch,phasedecoding,binarysearchwere introduced based on thelatesttheoreticalachievements ofthe algorithms.Finally,some applicationsbasedonSFTwereintroduced,anditsoutlookswerepresented. Keywords:SFT;spectrumpermutation;flatwindowfunction;subsampledFFT;Hashfunction 离散傅里叶变换(discreteFouriertransform, DFT)是信号处理中最基础的理论算法之一[1],计 算 DFT 最广泛的方法是1965年由 Cooley等[2]提 出的快速傅里叶变换,其时间复杂度为O(NlbN). 与直接计算 DFT 相比,FFT 大大简化了运算过程, 可将运算量缩短1~2个数量级,推动了信号处理技 术的革命性进展,被誉为二十世纪最具影响力的十 大算法之一[3].半个世纪以来,研究人员不断追求 FFT 算法的改进并提出了若干变体,如原地 FFT 算法和分裂基 FFT 算法[4]等,这些变体在某些硬件 实现上确有助益,但是复杂度并未得到明显改观. 人们开始将目光集中在信号自身特征上来,研究发 现,大量信号具有“稀疏性”:频域只有小部分大值 点,其余大部分点的值趋近于0.这一特征涵盖众多 应用领域:图像和语音信号[5]、机器学习理论[6]、布 尔函数 分 析[78]、压 缩 感 知[910]、宽 带 频 谱 感 知[11] 等.对于这种广泛存在且结构特殊的信号,能否实 现算法性能的突破呢? 首先进行这一研究的是1993年文献[12]中设 计的“哈达玛变换”快速算法,不久之后又出现了一 系列复数域上的研究[1316].这些算法在一定的稀疏 度范围内,性能超越 FFT,其时间复杂度大都是由

112 北京理工大学学报 第37卷 稀疏度K和1bN构成的多项式,如其中最快的是 文献[16]中的算法,时间复杂度为O(KbN)(指数 c>2).较大的指数项以及算法的复杂结构限制了 上述算法仅仅适应于“特别稀疏”的信号(K1),文献[17]对以上算法发展历程进行 图1分筐示意图 Fig. 1 Representation of the binning process 了较为详尽的描述.2012年MIT4位研究人员提 出一种更加简单有效的算法结构,在更广阔的适用 尽管SFT算法存在众多版本,但整体上都遵从 范围内性能超越FFT算法,重新掀起了研究的热图2所示的理论框架.分“筐”通过滤波器实现,为 潮,涌现出一大批算法理论成果[18-0、代码库[32了避免多个大值点出现在同一个“筐”中,首先对信 以及相关硬件设备[.本文的日的就是通过对这号频谱随机重排,使各大值点分散开.然后以带宽 些算法的研究,分析并阐明SFT算法原理,论述所为N/B的滤波器进行频段划分,频域的卷积使得每 涉及的主要技术难题,结合典型的算法将不同的重个频段内的点叠加在一起,再对频域以间隔N/B 构方法分门别类.另外,部分文献称这类算法为稀降采样,就能以极高的概率保证取样到所有大值点 疏快速傅里叶变换( sparse fast Fourier transfor 由B点FFT结果确定各大值点的值及所在“筐”的 SFFT)以凸显与FFT的联系,但这无关宏旨,后文位置.最后在保留的各“筐”中使用单频点重构的方 仍将沿用SFT的称法 法得到各大值点在原频谱中的坐标,继而求出较为 精确的近似值. 1理论框架及主要技术问题 有限列长为N的序列x(n)的DFT定义为 X(k)=∑x(n)W,k∈[0,N-1],(1) 图2SFT理论框架 式中WN=eN,假设信号x(n)是稀疏的,那么期 Fig 2 Theoretical framework of SFT 待一种“输出感知”的算法,能快速确定频域K个大1.1频谱随机重排 值点的坐标,得到频域近似值X(k),使其满足如下 对频谱重新排列最简单的方法是直接利用 2/2范数准则 DFT的两个基本性质 X-X‖2≤Cmin‖X-yl2,(2) ①位移性质:若x1(n)=x(n)Ww,则X1(k) 式中:参数C为误差系数;y为使|X-y2取最小X(k+b); 值的稀疏度为K的序列,即将X的K个最大值以 ②缩放性质:若x2(n)=x(an),则X2(k) 外的其它系数全部置0.该约束条件对估计谱X与X(ak) 信号实际谱X之间的误差做出了限定.另一个更 式中a-1表示关于模N的数论倒数,当a与N 加严格的约束条件是C。/2范数准则 互质,存在整数a1使得(-1)modN=1s.各参 ‖X-X|3≤X-y3/K+引‖x,(3)考文献中,频谱随机重排均通过以上两条性质实现 式中:精度参数c>0;6=1/N∞.2/2准则是对整以文献[22]为例 体估计误差的限定,而/2准则则保证频域每个系 p(n)=x(a(n-a))W,n∈[0,N-1], 数都能得到良好近似 SFT作为这样一种“输出感知”算法,其核心思式中:0与N互质;a、b为任意整数.对应的频域变 路是按照一定规则r(·)将信号频点投入到一组换关系为 筐”中(数量为B),如图1所示.因频域是稀疏的, P((k-b))=X(k)Ww,k∈[0,N-1]. 各大值点将依很高的概率在各自的“筐”中孤立存 (5) 在.将各“筐”中频点叠加,使N点长序列转换为B 可以看出,原频谱中k点通过变换后将出现在 点的短序列并作FFT运算,根据计算结果,忽略所(k-b)modN的位置上,根据数论中完全剩余系 有不含大值点的“筐”,最后根据对应分“筐”规则,设的概念,当σ与N互质,(k一b)modN依然是完备 计重构算法r(·)恢复出N点原始信号频谱 的,而Ww仅影响相位,故通过时域的变换实现了

稀疏度 K 和lbN 构成的多项式,如其中最快的是 文献[16]中的算法,时间复杂度为O(KlbcN)(指数 c>2).较大的指数项以及算法的复杂结构限制了 上述算法仅仅适应于“特别稀疏”的信号(K1),文献[17]对以上算法发展历程进行 了较为详尽的描述.2012年 MIT4位研究人员提 出一种更加简单有效的算法结构,在更广阔的适用 范围内性能超越 FFT 算法,重新掀起了研究的热 潮,涌现出一大批算法理论成果[1830]、代码库[3132] 以及相关硬件设备[3334].本文的目的就是通过对这 些算法的研究,分析并阐明 SFT 算法原理,论述所 涉及的主要技术难题,结合典型的算法将不同的重 构方法分门别类.另外,部分文献称这类算法为稀 疏快速傅里叶变换(sparsefastFouriertransform, SFFT)以凸显与 FFT 的联系,但这无关宏旨,后文 仍将沿用SFT 的称法. 1 理论框架及主要技术问题 有限列长为 N 的序列x(n)的 DFT 定义为 X(k)=∑ N-1 0 x(n)Wnk N ,k∈ [0,N -1],(1) 式中WN =e-j2π/N .假设信号x(n)是稀疏的,那么期 待一种“输出感知”的算法,能快速确定频域 K 个大 值点的坐标,得到频域近似值 X′ (k),使其满足如下 ℓ2/ℓ2范数准则 X -X′ 2 ≤C min K-sparsey X -y 2, (2) 式中:参数 C 为误差系数;y 为使 X-y 2 取最小 值的稀疏度为 K 的序列,即将 X 的K 个最大值以 外的其它系数全部置0.该约束条件对估计谱 X′ 与 信号实际谱X 之间的误差做出了限定.另一个更 加严格的约束条件是ℓ∞/ℓ2范数准则 X -X′ 2 ∞ ≤ε X -y 2 2/K +δ X 2 1,(3) 式中:精度参数ε>0;δ=1/NO(1).ℓ2/ℓ2准则是对整 体估计误差的限定,而ℓ∞/ℓ2准则则保证频域每个系 数都能得到良好近似. SFT 作为这样一种“输出感知”算法,其核心思 路是按照一定规则Γ(•)将信号频点投入到一组 “筐”中(数量为B),如图1所示.因频域是稀疏的, 各大值点将依很高的概率在各自的“筐”中孤立存 在.将各“筐”中频点叠加,使 N 点长序列转换为B 点的短序列并作 FFT 运算,根据计算结果,忽略所 有不含大值点的“筐”,最后根据对应分“筐”规则,设 计重构算法Γ-1(•)恢复出 N 点原始信号频谱. 图1 分筐示意图 Fig.1 Representationofthebinningprocess 尽管SFT 算法存在众多版本,但整体上都遵从 图2所示的理论框架.分“筐”通过滤波器实现,为 了避免多个大值点出现在同一个“筐”中,首先对信 号频谱随机重排,使各大值点分散开.然后以带宽 为N/B 的滤波器进行频段划分,频域的卷积使得每 个频段内的点叠加在一起.再对频域以间隔 N/B 降采样,就能以极高的概率保证取样到所有大值点, 由B 点 FFT 结果确定各大值点的值及所在“筐”的 位置.最后在保留的各“筐”中使用单频点重构的方 法得到各大值点在原频谱中的坐标,继而求出较为 精确的近似值. 图2 SFT理论框架 Fig.2 TheoreticalframeworkofSFT 1.1 频谱随机重排 对频 谱 重 新 排 列 最 简 单 的 方 法 是 直 接 利 用 DFT 的两个基本性质: ① 位移性质:若x1(n)=x(n)Wbn N ,则 X1(k)= X(k+b); ② 缩放性质:若 x2 (n)=x(σn),则 X2 (k)= X(σ-1k). 式中σ-1表示关于模 N 的数论倒数,当σ与 N 互质,存在整数σ-1使得(σσ-1)modN=1 [35].各参 考文献中,频谱随机重排均通过以上两条性质实现, 以文献[22]为例 p(n)=x(σ(n-a))Wσbn N ,n∈ [0,N -1], (4) 式中:σ与N 互质;a、b为任意整数.对应的频域变 换关系为 P(σ(k-b))=X(k)Wσak N ,k∈ [0,N -1]. (5) 可以看出,原频谱中k点通过变换后将出现在 σ(k-b)modN 的位置上,根据数论中完全剩余系 的概念,当σ与N 互质,σ(k-b)modN 依然是完备 的,而 Wσak N 仅影响相位,故通过时域的变换实现了 112 北 京 理 工 大 学 学 报 第 37 卷

第2期 仲顺安等:稀疏傅里叶变换理论及研究进展 113 对频谱的重新排序.其他变换形式如p(n)=器,其时域有效长度为O(B1bN6),参数B≥1, x(an+z)[21,均可由式(4)演化而来,不再赘述 1.2窗函数滤波器 δ>0,a>0,频域特性如下 由于每次重复或迭代循环中都涉及与滤波器相 ①G(k)=1,|k|≤(1-a)N/2B 乘,为了提高算法效率,要求该滤波器的时域和频域 ②G(k)=0,|k|≥N/2B, 都是稀疏的,即时域有效长度仅占小部分,频域除通 ③G(k)∈[0,1],(1-a)N/2B≤ 带以外均可以忽略. k|≤N/2B 文献[14,16]中使用的指示函数滤波器,时域 ④|G(k)-G(k)‖∞<δ 为矩形窗,对应频域为狄利克雷核序列,虽然也能实 在实现上,G(k)最理想的情况是矩形窗,根据 现分“筐”的功能,但是其阻带呈逆线性衰减,导致相对偶性,时域g(n)对应一个无限长的sinc序列,而 邻各筐的泄露较为严重 应用中的截断会造成边界点的不连续,引起频谱泄 文献[18]中提出一种脉冲序列滤波器,当露.一个有效的解决办法是,取两端平滑的窗函数 nmod n/p=0时,g(n)=1,n取其他值时,g(n) 与之时域相乘,如文献[21]中的高斯窗或文献[22] 0.它的优点是不存在泄露,频域是周期为p的脉冲中的切比雪夫窗,可以有效降低截断点的不连续性 序列但是要求参数p整除N,且算法需要多个不从而改善G(k)频域特性图3(a)所示为平坦窗函 同参数p,这对信号长度N提出了更多要求 数时域波形,有效长度仅占总长度的小部分,可降低 上述滤波器在很大程度上限制了SFT算法的所需样本量.图3(b)所示为平坦窗函数幅频特性, 发展,直到文献[21]中提出了一种平坦窗函数滤波通带内变得平坦,阻带得到指数衰减 2 a)平坦窗函数时域 (b)平坦窗函数频域 图3平坦窗函数滤波器设计 Fig 3 Design of the flat filtering windows function 13降采样FFT 对频域X(k)以相同间隔N/B降采样来获得所 y(n)=∑p(n+B),n∈[0,B-1 有的大值点 (7) Y(k)=P(kN/B),k∈[0,B-1].(6) 由此可见,经过信号时域的混叠,再做B点的 然而,无法直接对信号频谱进行操作,而是通过FFT运算就能得到包含所有大值点的频域降采样 时域的混叠来达到这一目的 序列.图4表示频域降采样的示意图,图4(a)中实 a)滤波后信号频域 (b)降采样后信号频域 图4频域降采样 Fig 4 Subsampling of spectrum

对频谱 的 重 新 排 序.其 他 变 换 形 式 如 p(n)= x(σn+τ)[21],均可由式(4)演化而来,不再赘述. 1.2 窗函数滤波器 由于每次重复或迭代循环中都涉及与滤波器相 乘,为了提高算法效率,要求该滤波器的时域和频域 都是稀疏的,即时域有效长度仅占小部分,频域除通 带以外均可以忽略. 文献[14,16]中使用的指示函数滤波器,时域 为矩形窗,对应频域为狄利克雷核序列,虽然也能实 现分“筐”的功能,但是其阻带呈逆线性衰减,导致相 邻各筐的泄露较为严重. 文献 [18]中 提 出 一 种 脉 冲 序 列 滤 波 器,当 nmodN/p=0时,g(n)=1,n取其他值时,g(n)= 0.它的优点是不存在泄露,频域是周期为p 的脉冲 序列,但是要求参数p 整除 N,且算法需要多个不 同参数p,这对信号长度 N 提出了更多要求. 上述滤波器在很大程度上限制了 SFT 算法的 发展,直到文献[21]中提出了一种平坦窗函数滤波 器,其时域有效长度为 O(B α lbN/δ),参数 B≥1, δ>0,a>0,频域特性如下: ① G′ (k)=1, k ≤(1-α)N/2B, ② G′ (k)=0, k ≥N/2B, ③ G′ (k)∈ [0,1], (1-α)N/2B ≤ k ≤N/2B, ④ G′ (k)-G(k) ∞ <δ. 在实现上,G(k)最理想的情况是矩形窗,根据 对偶性,时域g(n)对应一个无限长的sinc序列,而 应用中的截断会造成边界点的不连续,引起频谱泄 露.一个有效的解决办法是,取两端平滑的窗函数 与之时域相乘,如文献[21]中的高斯窗或文献[22] 中的切比雪夫窗,可以有效降低截断点的不连续性, 从而改善G(k)频域特性.图3(a)所示为平坦窗函 数时域波形,有效长度仅占总长度的小部分,可降低 所需样本量.图3(b)所示为平坦窗函数幅频特性, 通带内变得平坦,阻带得到指数衰减. 图3 平坦窗函数滤波器设计 Fig.3 Designoftheflatfilteringwindowsfunction 图4 频域降采样 Fig.4 Subsamplingofspectrum 1.3 降采样FFT 对频域X(k)以相同间隔N/B 降采样来获得所 有的大值点 Y(k)=P(kN/B), k∈ [0,B-1]. (6) 然而,无法直接对信号频谱进行操作,而是通过 时域的混叠来达到这一目的 y(n)= ∑ N/B-1 j=0 p(n+Bj), n∈ [0,B-1]. (7) 由此可见,经过信号时域的混叠,再做 B 点的 FFT 运算就能得到包含所有大值点的频域降采样 序列.图4表示频域降采样的示意图,图4(a)中实 第2期 仲顺安等:稀疏傅里叶变换理论及研究进展 113

114 北京理工大学学报 第37卷 线表示3个大值点经窗函数滤波器后的频谱,虚线O(bN)次循环中,对任意坐标k∈I=l1UI2U 表示降采样的位置,得到的结果如图4(b)所示 …∪I,若出现次数大于L/2,则将其归入集合Ⅰ 可以看出,第一个采样点在窗函数通带平坦区中,并认为集合Ⅰ包含所有目标频点坐标.对每 域,采样值与原谱线一致,这是最理想的情况,后两个k∈I,取L次循环得到的X(k)的中值作为最终 个频点采样位置分别出现在窗函数的上升沿和交叠的频率值,即 区域,采样值产生失真,应尽量避免.故要求随机重 X(k)= median({X(k)|r∈{1,2,…,L}}) 排使各大值点尽可能均匀分布,窗函数滤波器的上 哈希函数h。代表了分“筐”的规则,由频谱重排 升沿和下降沿尽量陡峭. 和窗函数共同决定.定位循环中通过h可将Z(k) 2重构方法 中保留的dK个大值点坐标反映射到原信号频谱, 得到dKN/B个坐标原象,最后联合L次循环结 通过频谱随机重排、窗函数滤波以及频域降采果,将其中出现频次最高的K个坐标值作为最终定 样将N点FFT转换为B点FFT,那么如何由B点位结果.对全部估值循环结果的实部和虚部分别取 FFT结果重构出信号完整的频谱呢?最典型的重中值作为最终大值点的频率值.据文献[21]的证 构方法有哈希映射法、混叠同余法、相位解码法和二明,该方法将以(1-1/N)的概率重构原信号频谱, 分查找法,下面将结合MIT算法最新理论成果,对满足式(2)范数准则,其时间复杂度为 这4种重构方法逐一进行介绍.首先约定信号长度 O(√ NKIb NlbN). N为2的整数次幂,如果不满足这一条件,可以通2.2混叠同余法 过人为添加若干0值点来达到 该方法最早由 Mansouri提出,适用于时域降 2.1哈希映射法 采样比较方便的信号,或者信号长度N为一系列互 该方法由 Hassanieh等在文献[21]中提出,得质整数的乘积,如N=70=2×5×7.其结构如下, 益于平坦窗函数滤波器的设计,相邻大值点之间的设参数M能整除N: 频谱泄露可以忽略,算法无需迭代,结构比较简单 ①随机选取偏移参数r∈[0,N/M-1] 为平衡分“筐”操作和短点FFT计算量,取参数B ②取z(n)=x(nN/M+r),n∈[0,M-1] √NK ③对z(n)作M点DFT运算得Z(k) ①频谱重排.取随机参数a,r(a为奇数)对信 ④将Z(k)的2K个最大值坐标归入集合T中 号进行重排p(n)=x(an+r),由式(5)知P(ahk) 时域的降采样对应频域的混叠,Z(k) X(k)wn: ②窗函数滤波.将p(n)通过1.2节所示滤波 X(k+nM)WNmM,由此式可知,X(k)中大 值点的坐标将以 kmod m的形式出现在集合T中 ③降采样.FFT将y(n)混叠 将参数M设置为多个不同值,得到一组线性同余方 程,通过中国剩余定理就可以求解出大值点的坐标 z(n)=∑y(n+B),n∈[o,B-1], 文献[21]中用此法对哈希映射法进行了优化,将 然后对z(m)作B点FFT运算,由式(7)知Z(k)=Y哈希映射的原象限定在集合modM∈T},有效降 (kN/B) 低反映射以及求中值的运算量,整个算法的时间复杂 ④哈希映射定义.哈希函数h。(k)= round度优化为O(√NK= lb nib n).弊端是无法避免碰撞 (aB/N), round表示4舍5人,定义偏移量: 问题,因为仅偏移参数随机取值,若j=j(modM), o,(k)=ak-h (k)(N/B) 即使经过频谱重排仍然有可=可(modM) ⑤定位循环.将Z(k)中dK个较大值的坐标k2.3相位解码法 归入集合J中,通过哈希反映射得到J的原像 该方法源于OFDM中的频偏估计方法,原 I={k∈[0,N-1]h2(k)∈J} 理是两个采样点之间的相位差与频点坐标成线性关 ⑥估值循环.对于每一个k∈I,计算对应的系.以单频信号为例,设x(n)=X(v)e2mN,∈ 频率值:X,(k)=Z(hn(k))W/G(on(k)) [0,N-1],则x(n+1)/x(n)=emN,通过相位 每一次定位循环得到一个坐标集合I,在L=2mxe/N就可以确定单频点坐标v.文献[22]使用

线表示3个大值点经窗函数滤波器后的频谱,虚线 表示降采样的位置,得到的结果如图4(b)所示. 可以看出,第一个采样点在窗函数通带平坦区 域,采样值与原谱线一致,这是最理想的情况,后两 个频点采样位置分别出现在窗函数的上升沿和交叠 区域,采样值产生失真,应尽量避免.故要求随机重 排使各大值点尽可能均匀分布,窗函数滤波器的上 升沿和下降沿尽量陡峭. 2 重构方法 通过频谱随机重排、窗函数滤波以及频域降采 样将 N 点 FFT 转换为B 点 FFT,那么如何由B 点 FFT 结果重构出信号完整的频谱呢? 最典型的重 构方法有哈希映射法、混叠同余法、相位解码法和二 分查找法,下面将结合 MIT 算法最新理论成果,对 这4种重构方法逐一进行介绍.首先约定信号长度 N 为2的整数次幂,如果不满足这一条件,可以通 过人为添加若干0值点来达到. 2.1 哈希映射法 该方法由 Hassanieh等在文献[21]中提出,得 益于平坦窗函数滤波器的设计,相邻大值点之间的 频谱泄露可以忽略,算法无需迭代,结构比较简单. 为平衡分“筐”操作和短点 FFT 计算量,取参数B≈ NK . ① 频谱重排.取随机参数σ,τ(σ为奇数)对信 号进行重排p(n)=x(σn+τ),由式(5)知 P(σk)= X(k)W-τk N ; ② 窗函数滤波.将p(n)通过1.2节所示滤波 器:y(n)=p(n)g(n); ③ 降采样.FFT 将y(n)混叠 z(n)= ∑ N/B-1 j=0 y(n+Bj),n∈ [0,B-1], 然后对z(n)作B 点 FFT 运算,由式(7)知Z(k)=Y (kN/B); ④ 哈 希 映 射 定 义.哈 希 函 数 hσ (k)=round (σkB/N),round表示4舍5入,定义偏移量: oσ(k)=σk-hσ(k)(N/B); ⑤ 定位循环.将Z(k)中dK 个较大值的坐标k 归入集合J 中,通过哈希反映射得到J 的原像: Ir={k∈[0,N-1]|hσ(k)∈J}; ⑥ 估值循环.对于每一个k∈Ir,计算对应的 频率值:Xr(k)=Z(hσ(k))Wτk N/G(oσ(k)). 每一次定位循环得到一个坐标集合Ir,在L= O(lbN)次 循 环 中,对 任 意 坐 标k∈I=I1 ∪I2 ∪ …∪IL,若出现次数大于 L/2,则将其归 入 集 合I′ 中,并认为集合I′ 包含所有目标频点坐标.对每一 个k∈I′ ,取L 次循环得到的X(k)的中值作为最终 的频率值,即: X(k)=median({Xr(k)|r∈ {1,2,…,L}}). 哈希函数hσ 代表了分“筐”的规则,由频谱重排 和窗函数共同决定.定位循环中通过hσ 可将Z(k) 中保留的dK 个大值点坐标反映射到原信号频谱, 得到dK N/B个 坐 标 原 象,最 后 联 合 L 次 循 环 结 果,将其中出现频次最高的 K 个坐标值作为最终定 位结果.对全部估值循环结果的实部和虚部分别取 中值作为最终大值点的频率值.据文献[21]的证 明,该方法将以 (1-1/N)的概率重构原信号频谱, 满足式(2)范数准则,其时间复杂度为 O( NKlbNlbN). 2.2 混叠同余法 该方法最早由 Mansour [13]提出,适用于时域降 采样比较方便的信号,或者信号长度 N 为一系列互 质整数的乘积,如 N =70=2×5×7.其结构如下, 设参数 M 能整除N: ① 随机选取偏移参数τ∈ [0,N/M -1]; ② 取z(n)=x(nN/M +τ),n∈ [0,M -1]; ③ 对z(n)作 M 点DFT 运算得Z(k); ④ 将Z(k)的2K 个最大值坐标归入集合T 中. 时 域 的 降 采 样 对 应 频 域 的 混 叠,Z(k)= ∑ N/M-1 n=0 X(k+nM)W-τ(k+nM) N ,由此式可知,X(k)中大 值点的坐标将以kmodM 的形式出现在集合T 中. 将参数 M 设置为多个不同值,得到一组线性同余方 程,通过中国剩余定理就可以求解出大值点的坐标. 文献[21]中用此法对哈希映射法进行了优化,将 哈希映射的原象限定在集合{kmodM ∈T},有效降 低反映射以及求中值的运算量,整个算法的时间复杂 度优化为O( 3NK2lbNlbN).弊端是无法避免碰撞 问题,因为仅偏移参数随机取值,若j=j ′ (modM), 即使经过频谱重排仍然有σj=σj ′ (modM). 2.3 相位解码法 该方法源于 OFDM 中的频偏估计方法[36],原 理是两个采样点之间的相位差与频点坐标成线性关 系.以单频信号为例,设x(n)=X(w)ej2πnw/N ,w∈ [0,N-1],则 x(n+1)/x(n)=ej2πw/N ,通 过 相 位 2πw/N 就可以确定单频点坐标w.文献[22]使用 114 北 京 理 工 大 学 学 报 第 37 卷

第2期 仲顺安等:稀疏傅里叶变换理论及研究进展 相位解码的方法对纯净的稀疏信号进行重构,此算 2/8-i1/2}, 构造信号:x(n)=x(2n)e-2×2m 对集合J中的每个元素 X(x)e2num,n∈[0,3],结合图5(b)观察: 1)令a=U(j)/U1(j); x(n+1)-ix(n)|<|x(n+1)+ir'(n) ‖)设不存在大值点碰撞,则a=e2mN,k= (12) a1( round( phase(a)N/2x)modN,其中 phase表 lx(n+1)-x(n)|<|x(n+1)+x(n) 示取相位; (13) ⅲ)与坐标k对应的频率值为U(k),Z(k) 假设式(12)成立且式(13)不成立,则对应( Z(k)+round(U(k)) 4)∈{1,2}属第二象限 对上述步骤进行O(bK)次迭代循环,将最终 继续构造信号:x(n)=x(2n)e-2m×1x2a/4 的Z(k)作为输入信号频谱X(k)的渐近估计值 (x)e2n(w-5)2,n∈[0,1],观察: 步骤②将N点FFT转换为B点FFT,相比于 x(n+1)-x(n)|<|x(n+1)+x(n) 从原信号中去除已知频点的影响,此处直接从FFT (14) 结果中将已知频点剔除,避免了IFT运算量.由 若上式成立,则-5=0,从而w=5. 式(4)可知,移位a点后进行同一重排变换其频域仅 存在相位差异Ws,步骤③利用a=0和a=1的两 次变换中相位差和频点坐标的线性关系,求得循环 中剩余大值点,并叠加至Z(k).据证明,经过 O(bK)次迭代,本算法以至少2/3的概率重构 X(k),时间复杂度为O(KlbN).尽管本算法较快, (b)N=4 但是鲁棒性极差,如x(n)=X(v)e2nmu/N+e(n),利 图5二分查找法重构单频信号图示 Fig 5 Recovering a single frequency via a binary search 用相位差计算α的前提条件是 e(n)|<|x(x)|/N 以上描述了二分查找法确定位置频点的过程 2.4二分查找法 该方法经过bN=3次循环,逐步缩小范围,直到判 对于普通含噪声信号,二分査找法具有更好的断出目标频点的位置.且二分查找法鲁棒性较好, 鲁棒性,结构相对比较复杂,本文将以最简单的单频即使信号受到小幅干扰,仍然可以判定出α=5.文 信号为例,阐述二分查找法的原理.设x(n)=献[22]采用二分查找法对普通含噪声信号进行重 X(ve)e2x/,t∈[o,7],是一个长度为8的单频信构,其时间复杂度为O( KIb NIb n/K) 号,频点坐标v未知,利用二分查找法确定频点位3SFT理论初步应用 置的方法如下 任取n∈[0,7],观察以下两式是否成立 自SFT提出以来,其应用潜力就初露端倪,然

相位解码的方法对纯净的稀疏信号进行重构,此算 法采用了迭代结构,初始化Z(k)=0,t代表循环次 数,结构如下: ① 随机参数σ为奇数,b∈[0,N-1],分别取 a=0和a=1,据式(4)作重排变换得p0 和p1; ② 对p0 和p1 分别执行以下运算,设输出分别 为U0 和U1: ⅰ)参数B=K/(2tβ); ⅱ)滤波后得y(n)=p(n)g(n),通过降采 样 FFT 计算Y(jN/B),j∈[0,B-1]; ⅲ)移除已知频点,计算Y′(jN/B)=Y(jN/ B)-(G′ (k)Z(σ-1b+k)Wa(σb+k) N )k=jN/B; ⅳ)返回Y′ (jN/B)的值. ③ 根据返回值U 求集合J={j:U(j)>1/2}, 对集合J 中的每个元素: ⅰ)令a=U0(j)/U1(j); ⅱ)设 不 存 在 大 值 点 碰 撞,则a=ej2πσk/N ,k= σ-1(round(phase(a)N/2π))modN,其中 phase表 示取相位; ⅲ)与坐标k 对应的频率值为U(k),Z(k)= Z(k)+round(U(k)). 对上述步骤进行 O(lbK)次迭代循环,将最终 的Z(k)作为输入信号频谱 X(k)的渐近估计值. 步骤②将 N 点 FFT 转换为B 点 FFT,相比于 从原信号中去除已知频点的影响,此处直接从 FFT 结果中将已知频点剔除,避免了IFFT 运算量.由 式(4)可知,移位a点后进行同一重排变换其频域仅 存在相位差异 Waσk N ,步骤③利用a=0和a=1的两 次变换中相位差和频点坐标的线性关系,求得循环 中剩 余 大 值 点,并 叠 加 至 Z(k).据 证 明,经 过 O(lbK)次 迭 代,本 算 法 以 至 少 2/3 的 概 率 重 构 X(k),时间复杂度为O(KlbN).尽管本算法较快, 但是鲁棒性极差,如x(n)=X(w)ej2πnw/N +ε(n),利 用 相 位 差 计 算 w 的 前 提 条 件 是 ε(n) < X(w)/N. 2.4 二分查找法 对于普通含噪声信号,二分查找法具有更好的 鲁棒性,结构相对比较复杂,本文将以最简单的单频 信号为例,阐述二分查找法的原理[37].设x(n)= X(w)ej2πnw/8,w∈[0,7],是一个长度为8的单频信 号,频点坐标 w 未知,利用二分查找法确定频点位 置的方法如下. 任取n∈[0,7],观察以下两式是否成立: ej2πw/8 -i < ej2πw/8 - (-i) , (8) ej2πw/8 -1 < ej2πw/8 - (-1) . (9) 结合图5(a)可知,式(8)对应于实轴上半区域, 式(9)对应于虚轴右半区域,且注意到x(n+1)= x(n)ej2πw/8,将上两式左右两边同时乘上 x(n) 得 x(n+1)-ix(n) < x(n+1)+ix(n) , (10) x(n+1)-x(n) < x(n+1)+x(n) . (11) 通过对以上两式的判断就可以确定所在象限, 假设以上两式均不成立,由式(10)可知,w∉{1,2, 3},由式(11)可知 w∉{0,1,7},故 对 应 第 三 象 限 w∈{4,5,6}. 构 造 信 号:x′ (n)= x (2n)e-j2π4×2n/8 = X(w)ej2π(w-4)n/4,n∈[0,3],结合图5(b)观察: x′ (n+1)-ix′ (n) < x′ (n+1)+ix′ (n) , (12) x′ (n+1)-x′ (n) < x′ (n+1)+x′ (n) . (13) 假设式(12)成立且式(13)不成立,则对应(w- 4)∈{1,2}属第二象限. 继续构 造 信 号:x″ (n)=x′ (2n)e-j2π×1×2n/4 = X(w)ej2π(w-5)n/2,n∈[0,1],观察: x″ (n+1)-x″ (n) < x″ (n+1)+x″ (n) . (14) 若上式成立,则w-5=0,从而w=5. 图5 二分查找法重构单频信号图示 Fig.5 Recoveringasinglefrequencyviaabinarysearch 以上描述了二分查找法确定位置频点的过程, 该方法经过lbN=3次循环,逐步缩小范围,直到判 断出目标频点的位置.且二分查找法鲁棒性较好, 即使信号受到小幅干扰,仍然可以判定出 w=5.文 献[22]采用二分查找法对普通含噪声信号进行重 构,其时间复杂度为O(KlbNlbN/K). 3 SFT理论初步应用 自SFT 提出以来,其应用潜力就初露端倪,然 第2期 仲顺安等:稀疏傅里叶变换理论及研究进展 115

116 北京理工大学学报 第37卷 而受限于复杂的算法结构和苛刻的前提条件,鲜有 ④算法理论推广,将SFT算法理论从傅里叶 较为出彩的应用.2012年以后,SFT技术取得长足域向其它变换域的推广,或者向更高维度的拓展,如 进展,为许多难题的解决带来了曙光.目前,SFT在文献[26,48]分别对二维和多维SFT进行了探索. 以下领域已有成功的应用 ①医学成像.大部分医学图像具有稀疏性,文5结束语 献[38]表明在2D相关振动光谱技术成像中利用 本文通过对SFT理论发展现状的研究,对其理 SFT算法能够有效减少扫描时间和截断伪影,有利论框架和典型算法结构进行了描述,介绍了所涉及 于减少对患者的危害并降低医疗成本 的关键技术,并对其中存在的问题进行了讨论,以期 ②频谱感知.是指通过信号检测和处理手段引起更多研究人员对这一领域的关注.SFT是一种 来获取无线网络中的频谱使用信息,对GHz带宽的较为新颖的理论,依然存在诸多问题有待研究,但它 检测会引起过高的采样率和超大的数据量,从而对是对传统傅里叶变换的一个极好的补充和完善,其 硬件设备造成巨大的挑战,文献[39]采用混叠同余研究成果可能对信号处理等领域产生重大影响」 类似方法,无需对全频带采样实现对稀疏信号的频 域进行计算,有效降低能耗和检测时间 参考文献: ③GPS同步,卫星传输时延较大,通常采用伪[1] Oppenheim A V. Schafer r w, Buck J R. Discrete-time 随机前导码来实现同步,文献[23]基于前导码较好 signal processing [m]. 2nd ed. Upper Saddle Rive 的自相关特性,IFFT结果仅含单根谱峰(K=1) NJ. USA: Prentice Hall, 1999 用SFT技术加速了帧同步处理过程 [2] Cooley J W, Tukey J W. An algorithm for the machine calculation of complex Fourier series[J].Mathematics ④其他应用.除以上应用外,SFT还被广泛应 of Computation, 1965, 19(90) 用到无线通信、光场摄影(、扩频捕获1、雷[37 Dongarra J, Sullivan F. Guest editors introduction to 达[等领域,同时在核磁共振(NMR)、光学相干层 the top 10 algorithms[J]. Computing in Science &.En- 析成像(OCT)、DNA测序、地震数据分析、无线天 gineering,2000,2(1):22-23 文学等方面显示出巨大应用潜力 [4] Rao K R, Kim D N, Hwang J.Fast Fourier transform-algorithms and applications [M]. Ist ed 4结论 Dordrecht: Springer Netherlands, 2010. 经过研究人员20年来的共同努力,SFT算法 [5] Anantha C, Gutnik V. Data driven signal processing: an 已经取得了一些显著成果,但是仍然存在很大提升 approach for energy efficient computing [c] l/ 空间.以下4个方面是当前研究的一些新趋势: Proceedings of International Symposium on Low Power ①进一步降低时间复杂度,尤其对含噪声信号 Electronics and Design. Monterey, CA, USA: IEEE 1996:347-352 来说,时间复杂度最低的是O( KIb NIb n/K),期待[6]ialN. Mansour y. Nisan N. Constant depth 将其优化为O(KbN),或在合理假设条件下证明 circuits,Fourier transform, and learnability [J] 其不存在 Association for Computing Machinery, 1993. 40(3) ②算法重构概率的提升,SFT算法只能以有限 607-620 的概率实现对信号频谱的重构,而简单通过增加循[7] Kahn J, Kalai G, Linial N. The influence of variables on 环或迭代次数来提高这一概率将牺牲时间复杂度 Boolean functions[C]//Proceedings of the 29th Annual 期待不以增加运算量为代价的前提下提高重构的 Symposium on Foundations of Computer Science. [S 概率 ③高速处理器或硬件平台下的实现,FT算法[8] O'Donnell R. Some topics in analysis of 已存在一系列快速代码库,如FFTw4 functions[C]//Proceedings of Annual ACM Symposium on Theory of Computing. Victoria, Canada: ACM CUFFT等,而目前SFT算法的源码3均基于 2008:569-578. x86CPUs,利用现代高速硬件平台或加速技术可使 [9] Donoho D L. Compressed sensing[J]. IEEE Transac- 性能得到进一步提升,如GPUs4, FPGAs[n, tions on Information Theory, 2006, 52(4): 1289-1306 DSPs以及其他多核处理器. [10] Candes E J, Romberg J, Tao T. Robust uncertainty

而受限于复杂的算法结构和苛刻的前提条件,鲜有 较为出彩的应用.2012年以后,SFT 技术取得长足 进展,为许多难题的解决带来了曙光.目前,SFT 在 以下领域已有成功的应用. ① 医学成像.大部分医学图像具有稀疏性,文 献[38]表明在 2D 相关振动光谱技术成像中利用 SFT 算法能够有效减少扫描时间和截断伪影,有利 于减少对患者的危害并降低医疗成本. ② 频谱感知.是指通过信号检测和处理手段 来获取无线网络中的频谱使用信息,对 GHz带宽的 检测会引起过高的采样率和超大的数据量,从而对 硬件设备造成巨大的挑战,文献[39]采用混叠同余 类似方法,无需对全频带采样实现对稀疏信号的频 域进行计算,有效降低能耗和检测时间. ③ GPS同步.卫星传输时延较大,通常采用伪 随机前导码来实现同步,文献[23]基于前导码较好 的自相关特性,IFFT 结果仅含单根谱峰(K=1),采 用SFT 技术加速了帧同步处理过程. ④ 其他应用.除以上应用外,SFT 还被广泛应 用到 无 线 通 信[40]、光 场 摄 影[41]、扩 频 捕 获[42]、雷 达[43]等领域,同时在核磁共振(NMR)、光学相干层 析成像(OCT)、DNA 测序、地震数据分析、无线天 文学等方面显示出巨大应用潜力. 4 结 论 经过研究人员20年来的共同努力,SFT 算法 已经取得了一些显著成果,但是仍然存在很大提升 空间.以下4个方面是当前研究的一些新趋势: ① 进一步降低时间复杂度,尤其对含噪声信号 来说,时间复杂度最低的是O(KlbNlbN/K),期待 将其优化为O(KlbN),或在合理假设条件下证明 其不存在. ② 算法重构概率的提升,SFT 算法只能以有限 的概率实现对信号频谱的重构,而简单通过增加循 环或迭代次数来提高这一概率将牺牲时间复杂度. 期待不以增加运算量为代价的前提下提高重构的 概率. ③ 高速处理器或硬件平台下的实现,FFT 算法 已 存 在 一 系 列 快 速 代 码 库,如 FFTW [44]、 CUFFT [45]等,而目前SFT 算法的源码[3132]均基于 x86CPUs,利用现代高速硬件平台或加速技术可使 性 能 得 到 进 一 步 提 升,如 GPUs [46],FPGAs [47], DSPs以及其他多核处理器. ④ 算法理论推广,将 SFT 算法理论从傅里叶 域向其它变换域的推广,或者向更高维度的拓展,如 文献[26,48]分别对二维和多维SFT 进行了探索. 5 结束语 本文通过对SFT 理论发展现状的研究,对其理 论框架和典型算法结构进行了描述,介绍了所涉及 的关键技术,并对其中存在的问题进行了讨论,以期 引起更多研究人员对这一领域的关注.SFT 是一种 较为新颖的理论,依然存在诸多问题有待研究,但它 是对传统傅里叶变换的一个极好的补充和完善,其 研究成果可能对信号处理等领域产生重大影响. 参考文献: [1]OppenheimAV,SchaferR W,BuckJR.Discrete-time signalprocessing[M].2nded.UpperSaddle River, NJ,USA:PrenticeHall,1999. [2]CooleyJW,TukeyJW.Analgorithmforthemachine calculationofcomplexFourierseries[J].Mathematics ofComputation,1965,19(90):297 301. [3]DongarraJ,SullivanF.Guesteditorsintroductionto thetop10algorithms[J].ComputinginScience& En￾gineering,2000,2(1):22 23. [4]Rao K R, Kim D N, Hwang J J. Fast Fourier transform-algorithmsand applications[M].1st ed. Dordrecht:SpringerNetherlands,2010. [5]AnanthaC,GutnikV.Datadrivensignalprocessing:an approach for energy efficient computing [C] ∥ ProceedingsofInternationalSymposiumonLowPower ElectronicsandDesign.Monterey,CA,USA:IEEE, 1996:347 352. [6]Linial N, Mansour Y, Nisan N. Constant depth circuits, Fourier transform, and learnability [J]. AssociationforComputing Machinery,1993,40(3): 607 620. [7]KahnJ,KalaiG,LinialN.Theinfluenceofvariableson Booleanfunctions[C]∥Proceedingsofthe29thAnnual Symposium onFoundationsofComputerScience.[S. l.]:IEEE,1988:68 80. [8]O’Donnell R. Some topics in analysis of boolean functions[C]∥ProceedingsofAnnualACMSymposium on Theory of Computing.Victoria,Canada:ACM, 2008:569 578. [9]DonohoD L.Compressedsensing[J].IEEE Transac￾tionsonInformationTheory,2006,52(4):1289 1306. [10]CandesEJ,RombergJ,TaoT.Robustuncertainty 116 北 京 理 工 大 学 学 报 第 37 卷

第2期 仲顺安等:稀疏傅里叶变换理论及研究进展 117 [22] incomplete frequency information [J]. IEEE Transac- ptimal sparse fourier transform[C]//Proceedings of tions on Information Theory, 2006, 52(2): 489-509 Annual ACM Symposium on Theory of Computing. [11] Lin M, Vinod A P. A new flexible filter bank for low [S..]:ACM,2012:563-577 complexity spectrum sensing in cognitive radios [J]. [23] Hassanieh H, Adid F, Katabi D, et al. Faster GPS via Journal of Signal Processing Systems, 2011, 62(2): the sparse fourier transform [c]// Proceedings of Annual Int Conf on Mobile Computing and [12] Kushilevitz E, Mansour Y. Learning decision trees Networking.[S.1.]:ACM,2012:353-364. using the Fourier spectrum[J]. SIAM Journal on Com- [24] Boufounos P, Cevher V. Gilbert A C, et al.What's puting,1993,22(6):1331-1348 the frequency, kenneth?: sublinear Fourier sampling [13] Mansour Y. Randomized interpolation and approxi off the grid [c]//Proceedings of Algorithms and tion of sparse polynomials [J]. SIAM Journal or Techniques Lecture Notes in Computer Science. omputing,l995,24(2):357-368 [S.1.]: Springer,2012:261 [14] Gilbert A C, Guha S, Indyk P, et al. Near-optimal [25] Lawlor D, Wang Y, Christlieb A. Adaptive sub-linear sparse Fourier representations via sampling [c] me Fourier algorithms [J]. Advances in Adaptive Proceedings of Annual ACM Symposium on Theory of Data Analysis, 2013, 5(1): 25 Computing. Montreal, Canada: ACM, 2002: 152- [26] Ghazi B, Hassanieh H, Indyk P, et al. Sample-optimal 161 average-case sparse Fourier transform in two [15] Akavia A, Goldwasser S, Safra S. Proving hard-core dimensions [C]//Proceedings of Communication predicates using list decoding [C]//Proceedings of Control, and Computing. [S.L. ] IEEE, 2013: 1258 IEEE Symposium on Foundations of Computer [27] Pawar S, Ramchandran K. Computing a k-sparse n- Science. Cambridge MA, USA: IEEE, 2003: 146 length discrete Fourier transform using at samples and O(klogk)complexity[C]// Proceedings of [16] Gilbert A C, Strauss M. Improved time bounds for IEEE International Symposium on Information near-optimal sparse Fourier representations [c] Theory.[S..]:IEEE,2013:464-468 Proceedings of International Society for Optical [28] Robin S, Saeid H, Martin V. A fast Hadamard Engineering. [S.I. ] SPIE,2005:59141-59115 transform for signals with sub-linear sparsity [c]// [17] Gilbert A C, Strauss M J, Tropp J A.A Proceedings of Communication, Control, and fast Fourier sampling: how to apply it to problems[J Computing. [S 1.]: IEEE, 2013:1250-I IEEE Signal Processing Magazine, 2008, 25(2): [29] Heider S, Kunis S, Munchen H Z, et al. A sparse Prony FFT[C]//Proceedings of the 10th Int Conf on [18]Wen ublinear-time ourier Sampling Theory and Applications. [S 1. ] Springer algorithms[J]. Foundations of Computational Mathe- 2013:572-575 matics,2010,10(3):303-338. [30] Indyk P, Kapralov M, Price E. Nearly) sample- [19] Akavia A. Deterministic sparse fourier approximation optimal sparse fourier transform[C]//Proceedings of via fooling arithmetic progressions[C]//Proceedings of Annual ACM-SIAM S Discrete he 23rd Conference on Learning Theory. [S 1]:Om- Algorithms. [S1]: ACM, 2014: 480-49 repress,2010:381-393. [31] Schumacher J. High Performance sparse fast fourier [20] Iwen M A. Improved approximation guarantees for transform[D]. Switzerland: ETH Zurich, 2013 sublinear-time Fourier algorithms [J]. Applied and [32] Wang C, Araya-Polo M, Chandrasekaran S, et al Computational Harmonic Analysis, 2013, 34(1) Parallel sparse FFT[C]//Proceedings ofthe Int Conf 57-82 for High Performance Computing, Networking [21] Hassanieh H, Indyk P, Katabi D, et al. Simple and Storage and Analysis. [S 1 ]: ACM,2013 practical algorithm for sparse fourier transform[C]/ [33] Abari O, Hamed E, Hassanieh H, et al. A0.75- Proceedings of Annual ACM-SIAM Symposium on million-point Fourier-transform chip for frequency Discrete Algorithms, [S. 1.]: ACM, 2012 sparse signals [C]//Proceedings of ISSCC. [S 1.] 1183-1194 IEEE,2014:458-460

principles:exact signal reconstruction from highly incompletefrequencyinformation[J].IEEE Transac￾tionsonInformationTheory,2006,52(2):489 509. [11]LinM,VinodAP.Anewflexiblefilterbankforlow complexityspectrum sensingincognitiveradios[J]. JournalofSignalProcessingSystems,2011,62(2): 205 215. [12]KushilevitzE, Mansour Y.Learning decisiontrees usingtheFourierspectrum[J].SIAMJournalonCom￾puting,1993,22(6):1331 1348. [13]Mansour Y.Randomizedinterpolation and approxi￾mationofsparsepolynomials[J].SIAM Journalon Computing,1995,24(2):357 368. [14]GilbertA C,GuhaS,IndykP,etal.Near-optimal sparse Fourierrepresentations via sampling[C]∥ ProceedingsofAnnualACM SymposiumonTheoryof Computing. Montreal,Canada: ACM,2002:152 161. [15]AkaviaA,GoldwasserS,SafraS.Provinghard-core predicatesusinglistdecoding[C]∥ Proceedings of IEEE Symposium on Foundations of Computer Science.Cambridge MA,USA:IEEE,2003:146 157. [16]GilbertA C,Strauss M.Improvedtimeboundsfor near-optimalsparse Fourier representations[C]∥ Proceedings of International Society for Optical Engineering.[S.l.]:SPIE,2005:59141 59115. [17]GilbertA C,StraussMJ,TroppJA.Atutorialon fastFouriersampling:howtoapplyittoproblems[J]. IEEE Signal Processing Magazine,2008,25 (2): 57 66. [18]Iwen M A. Combinatorial sublinear-time Fourier algorithms[J].FoundationsofComputationalMathe￾matics,2010,10(3):303 338. [19]AkaviaA.Deterministicsparsefourierapproximation viafoolingarithmeticprogressions[C]∥Proceedingsof the23rdConferenceonLearningTheory.[S.l.]:Om￾nipress,2010:381 393. [20]Iwen M A.Improvedapproximationguaranteesfor sublinear-time Fourieralgorithms[J].Applied and Computational Harmonic Analysis,2013,34 (1): 57 82. [21]HassaniehH,IndykP,KatabiD,etal.Simpleand practicalalgorithmforsparsefouriertransform[C]∥ Proceedings of Annual ACM-SIAM Symposium on Discrete Algorithms. [S.l.]: ACM, 2012: 1183 1194. [22]Hassanieh H,Indyk P, Katabi D,etal. Nearly optimalsparsefouriertransform[C]∥Proceedingsof AnnualACM Symposium on TheoryofComputing. [S.l.]:ACM,2012:563 577. [23]HassaniehH,AdidF,KatabiD,etal.FasterGPSvia thesparsefouriertransform [C]∥ Proceedings of Annual Int Conf on Mobile Computing and Networking.[S.l.]:ACM,2012:353 364. [24]BoufounosP,CevherV,GilbertA C,etal.What’s thefrequency,kenneth?:sublinearFouriersampling offthe grid[C]∥ Proceedings of Algorithms and Techniques Lecture Notes in Computer Science. [S.l.]:Springer,2012:261. [25]LawlorD,WangY,ChristliebA.Adaptivesub-linear timeFourieralgorithms[J].Advancesin Adaptive DataAnalysis,2013,5(1):25. [26]GhaziB,HassaniehH,IndykP,etal.Sample-optimal average-case sparse Fourier transform in two dimensions[C]∥ Proceedings of Communication, Control,andComputing.[S.l.]:IEEE,2013:1258. [27]PawarS,Ramchandran K.Computingak-sparsen￾lengthdiscrete Fouriertransform using at most4k samplesandO(klogk)complexity[C]∥Proceedingsof IEEE International Symposium on Information Theory.[S.l.]:IEEE,2013:464 468. [28]Robin S,Saeid H, Martin V. A fast Hadamard transformforsignals withsub-linearsparsity[C]∥ Proceedings of Communication, Control, and Computing.[S.l.]:IEEE,2013:1250 1257. [29]HeiderS,KunisS,München H Z,etal.Asparse PronyFFT[C]∥Proceedingsofthe10thIntConfon SamplingTheoryandApplications.[S.l.]:Springer, 2013:572 575. [30]IndykP,Kapralov M,Price E. (Nearly)sample￾optimalsparsefouriertransform[C]∥Proceedingsof Annual ACM-SIAM Symposium on Discrete Algorithms.[S.l.]:ACM,2014:480 499. [31]SchumacherJ.High Performancesparsefastfourier transform[D].Switzerland:ETHZurich,2013. [32]WangC,Araya-Polo M,Chandrasekaran S,etal. ParallelsparseFFT[C]∥ProceedingsoftheIntConf for High Performance Computing, Networking, StorageandAnalysis.[S.l.]:ACM,2013. [33]AbariO,Hamed E,Hassanieh H,etal.A 0.75- million-point Fourier-transform chip for frequency￾sparsesignals[C]∥ProceedingsofISSCC.[S.l.]: IEEE,2014:458 460. 第2期 仲顺安等:稀疏傅里叶变换理论及研究进展 117

118 北京理工大学学报 第37卷 [34] Yenduri P K, Rocca a Z, Rao A s,etal.Alow-[42]杜昌澔,韩航程,马瑛,等.基于稀疏傅里叶变换的深度 power compressive sampling time-based analog-to- 扩頻捕获与干扰抑制方法:中国,CN103427870A[P] digital converter [J]. IEEE Journal on Emerging and 2013-12-0 Selected Topics in Circuits and Systems, 2012, 2(3) Du Changhao, Han Hangcheng, Ma Ying, et al. A method for acquisition of spread spectrum an [35] Knuth D E. Art of computer programming, volume 2: interference suppression based on sparse FFT: China seminumerical algorithms [M]. [S. L ]:Addison- CN103427870A[P].2013-12-04.( in Chinese) Wesley, 1997 [43]陶然,刘升恒,张果,等,一种利用稀疏傅里叶变换计算 [36] Heiskala J, John Terry P D. OFDM wireless LANs:a 外辐射源雷达互模糊函数的方法:中国 theoretical and practical[M]. I ed. Indianapolis,IN CN103344948A[P].2013-10-09 USA: Sams Publishing, 2001 Rao, Liu Shengheng, Zhang Guo, et al. A fast al- [37] Gilbert A C, Indyk P, Iwen M, et al. Recent gorithms of cross-ambiguity function for passive radar developments in the sparse Fourier transform: a ased on external illuminato by using sparse FFT compressed Fourier transform for big data[J]. Signal hina,CN103344948A[P].2013-10-09.( in Chi- Processing Magazine, 2014, 31(5): 91-100 [38] Shi L, Andronesi O, Hassanieh H. MRS sparse-FFT: [44] Frigo M, Johnson S G. FFTW3 23[ EB/OL ]. reducing acquisition time and artifacts for in vivo 2d [2014-10-06].http://www.fftw.org/ orrelation spectroscopy [C]//Proceedi [45] Nvidia. CUDA CUFFT Library [EB/OL ] [201 International Society for Magnetic Resonance in 05-11].https://developer.nvidiacom/cufft Medicine Annual Meeting & Exhibition. [S. 1. ]: [46] Jiaxi H, Zhaosen W, Qiyuan Q, et al. Sparse fast ISMRM, 2013 Fourier transform on gpus and multi-core CPUs[C]// [39] Yenduri P K, Gilbert A C. Compressive, collaborative Proceedings of the 24th Int Symposium on Computer spectrum sensing for wideband cognitive radios[c]// Architecture and High Performance Computing. Proceedings of the 9th International Symposium on [S.I. ] IEEE Computer Society, 2012:83-91 Wireless Communication Systems. [S. L]: IEEE, [47] Agarwal A, Hassanieh H, Abari O, et al. High- throughput implementation of a million-point sparse [40] Hassanieh H, Lixin S, Bari O, et al. GHz-wide Fourier transform [C]//Proceedings of the 24th Int sensing and decoding using the sparse Fourier Conf on Field Programmable Logic and Applications. transform[c]//Proceedings of IEEE Conference [S.L.]:IEEE,2014:1-6 Computer Communications. [SL]: IEEE, 201 [48] Indyk P, Kapralov M. Sample-optimal Fourier 2256-2264. sampling in any constant dimension[C]//Proceedings [41] Shi L, Hassanieh H, Davis A, et al. Light field recon- of the 55nd Annual Symposium on Foundations of Computer Science. [S 1 ]:IEEE, 2014 domain [J]. ACM Transactions on Graphics,2014 (责任编辑:李兵) 31(1):1

[34]YenduriP K,Rocca A Z,Rao A S,etal.Alow￾power compressive sampling time-based analog-to￾digitalconverter[J].IEEEJournalonEmergingand SelectedTopicsinCircuitsandSystems,2012,2(3): 502 515. [35]KnuthDE.Artofcomputerprogramming,volume2: seminumericalalgorithms[M]. [S.l.]: Addison￾Wesley,1997. [36]HeiskalaJ,JohnTerryPD.OFDM wirelessLANs:a theoreticalandpractical[M].1ed.Indianapolis,IN, USA:SamsPublishing,2001. [37]Gilbert A C,Indyk P,Iwen M,et al. Recent developments in the sparse Fourier transform:a compressedFouriertransformforbigdata[J].Signal ProcessingMagazine,2014,31(5):91 100. [38]ShiL,AndronesiO,HassaniehH.MRSsparse-FFT: reducingacquisitiontimeandartifactsforinvivo2d correlation spectroscopy [C ] ∥ Proceedings of International Society for Magnetic Resonance in Medicine Annual Meeting & Exhibition. [S.l.]: ISMRM,2013. [39]YenduriPK,GilbertAC.Compressive,collaborative spectrumsensingforwidebandcognitiveradios[C]∥ Proceedingsofthe9thInternationalSymposium on Wireless Communication Systems. [S.l.]:IEEE, 2012:531 535. [40]Hassanieh H,Lixin S,Abari O,etal.GHz-wide sensing and decoding using the sparse Fourier transform[C]∥ProceedingsofIEEE Conferenceon Computer Communications. [S.l.]:IEEE,2014: 2256 2264. [41]ShiL,HassaniehH,DavisA,etal.Lightfieldrecon￾struction using sparsity in the continuous Fourier domain[J].ACM Transactionson Graphics,2014, 31(1):1. [42]杜昌澔,韩航程,马瑛,等.基于稀疏傅里叶变换的深度 扩频捕获与干扰抑制方法:中国,CN103427870A[P]. 2013 12 04. Du Changhao,Han Hangcheng,Ma Ying,etal.A method for acquisition of spread spectrum and interferencesuppressionbasedonsparseFFT :China, CN103427870A[P].2013 12 04.(inChinese) [43]陶然,刘升恒,张果,等.一种利用稀疏傅里叶变换计算 外 辐 射 源 雷 达 互 模 糊 函 数 的 方 法:中 国, CN103344948A[P].2013 10 09. TaoRao,LiuShengheng,ZhangGuo,etal.Afastal￾gorithmsofcross-ambiguityfunctionforpassiveradar basedon externalilluminato by usingsparse FFT: China,CN103344948A[P].2013 10 09.(inChi￾nese) [44]Frigo M, Johnson S G.FFTW3.23 [EB/OL]. [2014 10 06].http://www.fftw.org/. [45]Nvidia.CUDA CUFFT Library[EB/OL].[2014 05 11].https://developer.nvidia.com/cufft. [46]JiaxiH,Zhaosen W,Qiyuan Q,etal.Sparsefast Fouriertransformongpusandmulti-coreCPUs[C]∥ Proceedingsofthe24thIntSymposium onComputer Architecture and High Performance Computing. [S.l.]:IEEEComputerSociety,2012:83 91. [47]Agarwal A,Hassanieh H,Abari O,etal.High￾throughputimplementation ofa million-pointsparse Fouriertransform[C]∥Proceedingsofthe24thInt ConfonFieldProgrammableLogicand Applications. [S.l.]:IEEE,2014:1 6. [48]Indyk P, Kapralov M. Sample-optimal Fourier samplinginanyconstantdimension[C]∥Proceedings ofthe55nd AnnualSymposium on Foundations of ComputerScience.[S.l.]:IEEE,2014. (责任编辑:李兵) 118 北 京 理 工 大 学 学 报 第 37 卷