作者：桂。
時間：2018-01-06 14:00:25
鏈接：http://www.cnblogs.com/xingshansi/p/8214122.html
前言
對於數字接收來講，射頻域隨着帶寬的增加，AD、微波、FPGA資源的需求越來越高，但頻域開的越寬並不意味着頻譜越寬，有限信號內可認為信號在寬開頻域稀疏分布，最近較為流行的稀疏FFT（SFFT）是在傳統FFT的基礎上，利用了信號的稀疏特性，使得計算性能優於FFT。本文簡單記錄自己的理解。
一、稀疏FFT
主要是12年MIT的論文：Simple and Practical Algorithm for Sparse Fourier Transform。
核心思想：
SFT 作為這樣一種“輸出感知”算法,其核心思路是按照一定規則 Γ ( • )將信號頻點投入到一組“筐”中(數量為 B，通過濾波器實現 ) . 因頻域是稀疏的,各大值點將依很高的概率在各自的“筐”中孤立存在 . 將各“筐”中頻點疊加,使 N 點長序列轉換為 B點的短序列並作 FFT 運算,根據計算結果,忽略所有不含大值點的“筐”,最后根據對應分“筐”規則,設計重構算法 Γ -1 ( • )恢復出 N 點原始信號頻譜。
算法流程：
步驟一：選定一個sigma，實現數據頻域重排。
頻域重排主要是因為FFT之后，頻域大值集中在一起，需要將這一伙打散，保證頻域的稀疏性（sparse）。打散之后就能稀疏，依據定理（可見n越大，打散的可能性越大，從這里看，sigma與B還是有關系的，sigma體現了相鄰頻點的最小間隔，而B決定了每個bucket的寬度）：
這里是按一定概率（概率性）將頻譜重排，MIT重排思路：
對應代碼：
fs = 1000; f0 = 70; N = 256; t = [0:N-1]/fs; x = exp(1j*2*pi*t*f0) ; %%*************Step1:頻譜隨機重排************** sigma = 19; %inv(sigma) = 27 tao = 3; %permutation perm_num = mod([1:N]*sigma+tao,N)+1; pn = x(perm_num);
當然也有強制（確定性）將數據重排的思路，思路類似，只是實現方式不同：
sigma的選擇主要利用性質：
這里sigma^-1是sigma關於模N的逆元(數論倒數)。該點主要說明：變化前后的完備性（信號可重建）。
步驟二：加窗。
這里根據筐的多少（B），平分2pi頻域，即帶寬2pi/B，理想窗函數為矩形窗，但時域為無限長的sinc函數，需要加窗截斷，可選擇gausswin。
即win = sinc.*gausswin：
這里不局限於gausswin，滿足給定約束的窗函數均可：
步驟三：頻域抽取並作FFT
加窗之后，保證了頻譜不至於展寬嚴重，進一步保證了稀疏性。頻域降采樣：
等價於時域的混疊的傅里葉變換：
故直接在時域進行處理。處理完成后進行FFT運算。
該操作的理論基礎為：
不談加窗操作，x->p->y，可以看出x與y存在映射關系，而y的∑操作可以轉化為並行，這樣一來可以用多個低速率AD實現一個高速率AD的功能，前提是多個AD需要完全同步。
與上述降速率思想等價的是：如果各AD可做好同步，則多個現有AD的能力，可以做出現有AD難以完成的事。
步驟四：哈希映射
哈希映射的線索為：最終觀察的數據Y(k)【即y(n)】->P(k)【即p(n)】->X(k)【即x(n)】：
第一步（y -> p）轉化的理論依據：
第二步（p -> x）轉化的理論依據：
序列重排后的頻域變換為：
這兩步解釋了算法step4的參數定義：
但實際中並非第一步提到的理想情況，實際情況是：
因此存在一個頻譜重建的成功概率問題：
4章節的后半部分主要在證明這個概率問題。
步驟五~六的主要目的，引用原文的說法：
Here are two versions of the inner loop: location loops and estimation loops. Location loops are given a parameter d, and output a set I ⊂ [n] of dkn/B coordinates that contains each large coefficient with“good” probability. Estimation loops are given a set I ⊂ [n] and estimate x I such that each coordinate is estimated well with “good” probability.
步驟五：定位循環
d是一個參數，存在約束
原文取值為：
該步驟的主要操作為：
將 Z ( k )中 dK 個較大值（從大到小依次排列）的坐標（ k）歸入集合 J 中,通過哈希反映射得到 J 的原像：
這樣便得到包含原始頻譜坐標的集合，迭代L次：
迭代的目的？
步驟六：估值循環
對於每個k∈Ir，計算頻率信息：
步驟5~6主要操作：
至此完成了稀疏FFT（sparse FFT）的整個運算過程，記作sfft_v01。
改進版都是依次大框架進行，不同點主要有三個方面：1）重排的實現思路；2）頻譜的重建思路。不再一一展開。
原文示例：
該示例給出了步驟5~6的直觀解釋：
二、問題記錄
仿真驗證：n = 256點頻信號，sigma = 19，則inv_sigma = 27，假設生成點頻信號：
圖1為重排的信號，圖2為重排加窗(此處窗不夠理想，看到長尾，論文中強調加窗，作用在於把尾巴剁掉。)，圖3為原始信號頻譜，圖4為降采樣之后的頻譜。可以看出重排容易引入諧波：以sigma為周期。
原始頻譜位置為k=10-1 = 9（下標從0開始），重排后對應頻譜位置為mod(sigma*k,n) = 171 = 172-1，與理論分析一致。即原始k點重排后位置：
即時域信號，x(n)，重排序號perm_num = mod([1:n]*sigma,n)+1，則對應重排后的傅里葉變換，頻譜順序與原信號頻譜X，經過重排序號perm_num1 = mod([1:n]*inv_sigma,n)+1，位置始終相差1（數論倒數的性質決定）:mod(sigma*inv_sigma,n)=1