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

第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).较大的指数项以及算法的复杂结构限制了 上述算法仅仅适应于“特别稀疏”的信号(K<6(N/ lbN),a>1),文献[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).较大的指数项以及算法的复杂结构限制了 上述算法仅仅适应于“特别稀疏”的信号(K<Θ(N/ lbaN),a>1),文献[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-i<|e2m-(-i)|,(8) 法采用了迭代结构,初始化Z(k)=0,t代表循环次 em6-1|<|e2m/-(-1).(9) 数,结构如下: 结合图5(a)可知,式(8)对应于实轴上半区域 ①随机参数a为奇数,b∈[0,N-1],分别取式(9)对应于虚轴右半区域,且注意到x(n+1) 0和a=1,据式(4)作重排变换得p和p1; x(n)e2m/8,将上两式左右两边同时乘上|x(n)得 ②对p和p1分别执行以下运算,设输出分别 x(n+1)-ix(n)<|x(n+1)+i(n), 为U。和U1 (10) 参数B=K/(2B); lx(n+1)-x(n)|<|x(n+1)+x(n) ‖)滤波后得y(n)=p(n)g(n),通过降采样 (11) FT计算Y(N/B),j∈[0,B-1] 通过对以上两式的判断就可以确定所在象限, ⅲ)移除已知频点,计算Y(jN/B)=Y(jN/假设以上两式均不成立,由式(10)可知,∈{1 B)-(G(k)Z(o-1b+k)wo+)Ie-iN/B 3},由式(11)可知v∈{0,1,7},故对应第三象限 )返回YGjN/B)的值 ∈{4,5,6} ③根据返回值U求集合J=5:|U(j)|>1/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