NUFFT的計算難點在哪裡?

01-11

NUFFT （NonUniform Fast Fourier Transform）的計算難是相對於 FFT （Fast Fourier Transform）而言的。NUFFT和FFT屬於 DFT（Discrete Fourier Transform）處理不同離散信號的快速計算演算法。FFT是處理等間隔採樣的離散信號，而NUFFT是處理採樣間隔不固定（採樣密度變化）的離散信號。

前難點. NUFFT 的計算比 FFT慢。NUFFT包括兩個步驟：卷積插值（gridding）和FFT，其計算複雜度由卷積插值步驟決定。最近幾年對卷積插值快速計算的研究，有包括優化卷積窗函數（在保證插值精度的前提下縮小窗寬，降低過採樣比等），有基於openmp加速的NFFT庫和基於GPU加速的gpuNUFFT庫等。所以計算慢已經不是NUFFT的突出問題。

現難點. NUFFT計算結果的理解，可逆性，和逆變換的計算問題。

先看FFT變換：在滿足奈奎斯特採樣定理條件下，其對一個帶寬有限的一維連續信號做等間隔採樣得到長度為N的離散信號x[n]，其離散傅里葉變換記作矩陣向量相乘的形式為：

f = F*x,

其中F為大小為N x N 的傅里葉變換矩陣，f[k] 是該信號的離散頻譜，在K個等間隔頻率點上，一般K = N。該變換有FFT演算法可以快速計算（FFT的優點1）。從這個頻譜再經過inverse FFT變換：

x = F^(-1) * f

仍能夠得到原信號（FFT的優點2），F^(-1)表示F的逆。在數學上，F矩陣是可逆的，並且有F矩陣的共軛轉置等於F矩陣的逆，也即F^(-1) = F^(H) ，而F^(H) 同樣可以用FFT演算法快速計算（FFT的優點3）.

所以FFT變換是可逆的，並且正逆變換都有快速計算演算法。

但是對同一個連續信號做非等間隔採樣（變密度採樣）得到長度仍為N的離散信號t[n]，做NUFFT變換並假設得到的結果頻譜y的頻率點是等間隔的，其矩陣向量乘的形式為：

y = E * t

那麼這個頻譜 y 和 f 是不等價的（問題1）。從這個頻譜 y 經過逆變換（t1 = E^(-1) * y）得到的 t1等不等於t，以及這個逆變換E^(-1) 與 E^(H) (E矩陣的共軛轉置)之間的關係是什麼都成了問題。

變換結果的理解問題：目前仍屬於一個研究問題。

變換結果的可逆性問題：可逆的存在性問題可以通過計算E矩陣的秩（rank）獲得，但該方法限於處理較小的E矩陣。逆矩陣的求解目前有迭代和非迭代法。非迭代法：引入採樣密度補償函數W，使得 E*W*t 的結果近似等於 F*x，或者使得 E^(H) * W * y 近似等於t。W的計算有多種方法，但是近似精度一直是個問題。迭代法：使用共軛梯度演算法迭代求解 E矩陣的偽逆，計算量是迭代次數乘以2，再乘以一次NUFFT的計算量。

所以NUFFT的現難點在於變換的結果怎麼理解，變換結果是否可逆，可逆的條件，NUFFT變換及逆變換的快速計算方法。

總結：NUFFT的難已不在計算，而是難在NUFFT的計算結果怎麼理解和怎麼用。

例子：磁共振成像里使用Radial或Spiral軌跡採樣，就是變密度採樣和NUFFT演算法的一個應用。在磁共振中，Radial採樣有自己的優勢。但是如何求逆重建圖像，重建圖像的可用性，偽影的產生及抑制，如何設計合適的採樣模式（採樣密度的分布情況）等都是上述NUFFT演算法難點的具體化。

由於uniform Fourier transform有快速演算法，可以將計算量從N^2減少到N log（N），就是俗稱的FFT，而非格點分布的信號則不能使用FFT，所以就產生了對應的快速演算法，名為NUFFT（non uniform fast Fourier transform）。其實思想很簡單，就是先把非格點分布的信號「放」在格點上，然後再使用FFT。第二步沒啥好說的，難點在第一步。

第一步可以有很多種方法來實現。目前通用的方法由以下幾步：

1. density compensation

2. interpolate onto a finer grid using some core function

3. Deapodization

1. 由於非均勻分布的信號，如果直接將它們「放」在附近的格點上，那麼會造成格點上信號分布的不均勻，所以需要一個密度補償。

2. 通常的做法是差值到一個更為密集的格點上，然後在進行完FFT操作之後再crop到需要的大小，這樣做可以避免aliasing。之前通用的做法是2倍的grid，但後來有文章指出其實1.1～1.3倍的就夠了。具體數字記不太清了。core function和density conpensation主導了NUFFT的精度。完美core是一個無限長的sinc function，但不實際，所以通用的做法是用一個短一點的方程。目前很好的一個近似core是Kaiser-Bessel。

3. 由於在第二步使用了一個core，所以在這裡要補償一個core的Fourier transform，不然密度會出問題。

大概就是這些，寫到一半不想寫了，以後有閑了再補充。。。