稀疏傅里葉變換原理說明（一）

前言

因為某些情況，鑽研了幾個月的稀疏傅里葉變換，總算從“一無所知”到“對其有了一些了解”，因此寫下這篇博文以便參考。一是希望給剛開始研究稀疏傅里葉變換的網友一些經驗，少走彎路；二是由於博主並非通信專業出身，對數字信號處理幾乎一竅不通，導致即使我知道了稀疏傅里葉變換的大致原理，仍然無法徹底搞懂這個算法，程序中的各種細節也看不明白，因此希望和廣大網友相互交流借鑒。

由於水平有限，博文之中可能存在解釋不當的地方甚至是錯誤，歡迎各位指正！

參考代碼鏈接：bernielampe1/sparse_fft: Implementation of sparse fast Fourier transform (sFFT). (github.com)

一、稀疏傅里葉變換簡介

稀疏傅里葉變換，英文名為Sparse Fourier Transform，一般縮寫為SFT或者SFFT（Sparse Fast Fourier Transform），下文中都稱之為SFT。SFT是一種利用信號自身特點得到信號頻譜的算法，是對信號處理中最經典的算法之一——離散傅里葉變換（縮寫為DFT，下同）的一種改進。所謂“信號自身特點”，對一維信號而言，一般指的是其頻率范圍內只有少部分頻率分量幅值較大，其他大部分頻率分量的幅值都趨近或者等於0。舉個例子：一個信號在0~1000Hz的范圍內，只有50Hz和120Hz這兩個頻率上具有能量，其他頻率上幾乎沒有能量。這一特點也被稱為“稀疏性”。

注：一個只有在50Hz和120Hz上有“有用信號”的稀疏信號

稀疏性一般使用稀疏度來表示。稀疏度（用K表示）即有用頻率成分的個數。以上圖為例，稀疏度K=2。

SFT最先是由MIT在2012年的論文《Simple and Practical Algorithm for Sparse Fourier Transform》中提出來的，該博客也主要是基於該論文進行的研究。

二、稀疏傅里葉變換流程

SFT的理論框架如圖。

這里對各個部分進行簡要說明：

• 頻譜重排：利用DFT中頻域和時域的性質，對時域序列進行操作，使得信號的頻譜重新排列
• 濾波：使用提前設計好的濾波器，將之與重排后的信號相乘，起到擴展譜線、截短序列的作用
• 降采樣FFT：仍然是利用DFT的性質，通過對時域序列求和，實現對頻譜的采樣
• 定位：利用哈希逆變換，從降采樣后得到的頻譜反推出原信號譜線所在的位置
• 求值：利用濾波器的頻譜和估計得到的位置，反推出原信號譜線的幅值

下面會對各個部分進行詳細解釋。

三、流程中各部分詳解

3.1 頻譜重排

頻譜重排主要是利用了DFT的性質：

1. 位移性質：若時域上有：$$ x_{1}=x(n)W_{N}^{bn}$$ 則頻域上有：$$ X_{1}(k)=X(k+b) $$
2. 縮放性質：若$$ x_{2}=x(\sigma n)$$ 則頻域上有：$$ X_{2}(k)=X(\sigma^{-1} k) $$

其中σ^-1並非和σ互為倒數，而是一個整數，該整數滿足關系：$$\sigma \sigma^{-1} mod(n) = 1$$

也就是σ^-1與σ的乘積對n取模，所得余數等於1。此時稱σ^-1是σ關於模n下的逆元（博客中n和N含義相同，n論文中常用，N個人更喜歡用）。舉兩個例子：

7在模24下的逆元就是本身，因為$$7\times 7\div 24=2\cdots 1$$

15在模128下的逆元是111，因為$$15\times 111\div 128=13\cdots 1$$

當然並非所有數都存在逆元，不過這不在討論范圍之內。

縮放性質是SFT中重要的理論基礎，其可以使用matlab進行驗證。這里需要注意的是，超出了序列范圍之外的下標，需要對序列長度求模進行約束。

N = 128;
F_signal = zeros(N, 1);
F_signal(60) = 5;
signal = ifft(F_signal); % 利用反變換得到信號序列

p_signal = zeros(N, 1);

% 進行頻譜重排，其中重排因子為 15，其關於模 128 的倒數為 111
for i = 0:N-1
    p_signal(i+1) = signal(mod(i*15, N)+1);
end

F_p_signal = fft(p_signal); % 對新序列進行傅里葉變換，查看其頻譜

subplot(2,1,1);
stem(F_signal);
subplot(2,1,2);
stem(abs(F_p_signal));

所得結果如圖：

由於matlab下標不能從0開始，因此上半部分的譜線實質對應着序號59。以15作為重排因子，128作為模，所得新信號的譜線對應位置為117，滿足：

$$\left\{\begin{matrix}
59\times 15\div 128=6\cdots 117
\\ 
117\times 111\div 128=101\cdots 59
\end{matrix}\right.$$

3.2 濾波

雖然叫濾波，但這個的濾波過程並非是傳統意義上的濾出信號中的有用頻率成分。首先需要知道的是，在時域上兩個信號進行卷積操作，等於在頻域上對這兩個信號進行相乘操作。反過來，在時域上對兩個信號的對應元素進行相乘，則對應於在頻域上對這兩個信號進行卷積。而這里的濾波，指的就是后一種方式。下圖就顯示了這兩種方式的區別。

N = 256;
omega = 50;

% 生成信號和幅值
F_signal = zeros(N, 1);
F_signal(117) = 4;
F_signal(210) = 1;
signal = ifft(F_signal);

% 生成矩形窗濾波器
F_box_car = zeros(N, 1);
F_box_car((N-omega)/2 : (N+omega)/2) = 1;
box_car = ifft(F_box_car);
F_box_car = fft(box_car);

% 相乘
times_b_s = signal .* box_car;
F_times_b_s = fft(times_b_s);

% 卷積
conv_b_s = conv(signal, box_car, 'same');
F_conv_b_s = fft(conv_b_s);

% 繪圖
subplot(4,1,1);
plot(abs(F_signal)); title('原信號的頻譜');
subplot(4,1,2);
plot(abs(F_box_car)); title('濾波器頻譜');
subplot(4,1,3);
plot(abs(F_times_b_s)); title('時域上相乘后信號的頻譜');
subplot(4,1,4);
plot(abs(F_conv_b_s)); title('時域上卷積后信號的頻譜');

 

第四個子圖中展示的行為就是我們常見的濾波器的操作。在第三個子圖中，原信號的頻率分量被濾波器展開並移動，而從幅值的大小可以大致猜出發生了怎樣的變化。

使用時域相乘這種方法的要點有兩個：

1. 同一頻率分量經濾波器擴展后，幅值相同（假如濾波器理想），因此在頻域上采樣得到的幅值可以代表原信號頻率成分的幅值；
2. 降低了運算量，比如原信號序列長為4096，但濾波器序列可以僅為256，只需取出信號序列中的256個點與濾波器進行相乘即可。

值得一提的是，濾波器的性能直接決定了SFT的可靠性，而濾波器問題也是困擾我的問題之一。理想情況下，直接使用矩形窗濾波器就能完成所有事情，但理想的矩形窗濾波器對應着一個時域上無限長的sinc序列，這顯然無法實現。現實中，我們常用過渡帶平滑的濾波器與理想矩形窗相乘或卷積來改善矩形窗性能，使之在一定序列長度下，幅頻特性盡量接近理想矩形窗的幅頻特性。處理后的矩形窗被稱為“平坦窗”。

如何實現平坦窗的自動化生成，一直是我在考慮的問題。論文中關於這一段的描述如下：

Let B be a parameter that divides n, to be determined later. Let G be a (1/B, 1/(2B), δ, ω) flat window function for some δ and ω= O(Blog(n/δ)). We will have δ≈1/n^c, so one can think of it as negligibly small.

B是分筐數，后面會提。B一般滿足關系

$$B\approx \sqrt{nk}$$

以序列長度n=8192、稀疏度K=7為例，B=239，由於B應該能整除n，因此取與之最接近的256。按照上面的公式，取δ≈1/n^c≈1/n，則平坦窗寬度ω=256*log(8192²)≈4614，顯然與時域相乘中的要點（2）相沖突。即使將ω進行縮小，縮小到何種程度也是一個問題。另一方面，平坦窗由兩個窗卷積得到，一般選擇矩形窗和高斯窗。如何利用已知參數（分筐數B、稀疏度K、序列總點數N等）計算出合適的參數來生成矩形窗和高斯窗，兩者的卷積能滿足(1/B, 1/(2B), δ, ω)這一條件，也是一個非常麻煩的問題。

---

 

由於篇幅限制，后續內容放在了“稀疏傅里葉變換原理說明（二）”和“稀疏傅里葉變換原理說明（三）”當中。