作者：桂。

時間：2018-01-06??14:00:25

鏈接：http://www.cnblogs.com/xingshansi/p/8214122.html?

---

前言

　　對于數(shù)字接收來講，射頻域隨著帶寬的增加，AD、微波、FPGA資源的需求越來越高，但頻域開的越寬并不意味著頻譜越寬，有限信號內(nèi)可認為信號在寬開頻域稀疏分布，最近較為流行的稀疏FFT（SFFT）是在傳統(tǒng)FFT的基礎上，利用了信號的稀疏特性，使得計算性能優(yōu)于FFT。本文簡單記錄自己的理解。

源代碼

一、稀疏FFT

主要是12年MIT的論文：Simple and Practical Algorithm for Sparse Fourier Transform。

核心思想：

SFT 作為這樣一種“輸出感知”算法,其核心思路是按照一定規(guī)則 Γ ( ? )將信號頻點投入到一組“筐”中(數(shù)量為 B，通過濾波器實現(xiàn) ) . 因頻域是稀疏的,各大值點將依很高的概率在各自的“筐”中孤立存在 . 將各“筐”中頻點疊加,使 N 點長序列轉(zhuǎn)換為 B點的短序列并作 FFT 運算,根據(jù)計算結(jié)果,忽略所有不含大值點的“筐”,最后根據(jù)對應分“筐”規(guī)則,設計重構算法 Γ ^-1 ( ? )恢復出 N 點原始信號頻譜。

算法流程：

步驟一：選定一個sigma，實現(xiàn)數(shù)據(jù)頻域重排。

頻域重排主要是因為FFT之后，頻域大值集中在一起，需要將這一伙打散，保證頻域的稀疏性（sparse）。打散之后就能稀疏，依據(jù)定理（可見n越大，打散的可能性越大，從這里看，sigma與B還是有關系的，sigma體現(xiàn)了相鄰頻點的最小間隔，而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);

　　當然也有強制（確定性）將數(shù)據(jù)重排的思路，思路類似，只是實現(xiàn)方式不同：

sigma的選擇主要利用性質(zhì)：

這里sigma^-1是sigma關于模N的逆元(數(shù)論倒數(shù))。該點主要說明：變化前后的完備性（信號可重建）。

?步驟二：加窗。

?這里根據(jù)筐的多少（B），平分2pi頻域，即帶寬2pi/B，理想窗函數(shù)為矩形窗，但時域為無限長的sinc函數(shù)，需要加窗截斷，可選擇gausswin。

即win = sinc.*gausswin：

這里不局限于gausswin，滿足給定約束的窗函數(shù)均可：

步驟三：頻域抽取并作FFT

加窗之后，保證了頻譜不至于展寬嚴重，進一步保證了稀疏性。頻域降采樣：

等價于時域的混疊的傅里葉變換：

故直接在時域進行處理。處理完成后進行FFT運算。

該操作的理論基礎為：

不談加窗操作，x->p->y，可以看出x與y存在映射關系，而y的∑操作可以轉(zhuǎn)化為并行，這樣一來可以用多個低速率AD實現(xiàn)一個高速率AD的功能，前提是多個AD需要完全同步。

與上述降速率思想等價的是：如果各AD可做好同步，則多個現(xiàn)有AD的能力，可以做出現(xiàn)有AD難以完成的事。

步驟四：哈希映射

哈希映射的線索為：最終觀察的數(shù)據(jù)Y(k)【即y(n)】->P(k)【即p(n)】->X(k)【即x(n)】：

第一步（y -> p）轉(zhuǎn)化的理論依據(jù)：

第二步（p -> x）轉(zhuǎn)化的理論依據(jù)：

序列重排后的頻域變換為：

這兩步解釋了算法step4的參數(shù)定義：

但實際中并非第一步提到的理想情況，實際情況是：

因此存在一個頻譜重建的成功概率問題：

?4章節(jié)的后半部分主要在證明這個概率問題。

步驟五~六的主要目的，引用原文的說法：

　　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.

?步驟五：定位循環(huán)

?d是一個參數(shù)，存在約束

原文取值為：

?該步驟的主要操作為：

將 Z ( k )中 dK 個較大值（從大到小依次排列）的坐標（ k）歸入集合 J 中,通過哈希反映射得到 J 的原像：

這樣便得到包含原始頻譜坐標的集合，迭代L次：

迭代的目的？

?

步驟六：估值循環(huán)

?對于每個k∈Ir，計算頻率信息：

步驟5~6主要操作：

至此完成了稀疏FFT（sparse FFT）的整個運算過程，記作sfft_v01。

改進版都是依次大框架進行，不同點主要有三個方面：1）重排的實現(xiàn)思路；2）頻譜的重建思路。不再一一展開。?

原文示例：

該示例給出了步驟5~6的直觀解釋：

?

二、問題記錄

仿真驗證：n = 256點頻信號，sigma = 19，則inv_sigma = 27，假設生成點頻信號：

圖1為重排的信號，圖2為重排加窗(此處窗不夠理想，看到長尾，論文中強調(diào)加窗，作用在于把尾巴剁掉。)，圖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，經(jīng)過重排序號perm_num1 =??mod([1:n]*inv_sigma,n)+1，位置始終相差1（數(shù)論倒數(shù)的性質(zhì)決定）:mod(sigma*inv_sigma,n)=1

總結(jié)

以上是生活随笔為你收集整理的稀疏傅里叶变换（sparse FFT）的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。