这篇文档主要介绍bluestein算法(又叫Chirp-Z变换)的实现原理以及GPU上实现的思路。

算法原理

离散信号的傅里叶变换可以写成以下公式:

$$ X_k=\sum_{n=0}^{N-1}x(n)e^{-\frac{2i\pi}{N}nk} $$

观察这个表达式,可以看到$nk$这一项可以用平方式展开:

$$ (n-k)^2=n^2-2nk+k^2 => nk=-\frac{(n-k)^2-n^2-k^2}{2} $$

将$nk$这一项带入原来的公式:

$$ e^{-\frac{2i\pi}{N}nk}=e^{\frac{2i\pi}{N}\frac{(n-k)^2-n^2-k^2}{2}}=e^{\frac{\pi i(n-k)^2}{N}}e^{-\frac{\pi in^2}{N}}e^{-\frac{\pi ik^2}{N}} $$

因此,原来的傅里叶变换公式可以写为:

$$ X_k=e^{-\frac{\pi i k^2}{N}}\sum_{n=0}^{N-1}{x(n)}e^{\frac{\pi i(n-k)^2}{N}}e^{-\frac{\pi in^2}{N}} $$

设$A_n=x(n)e^{\frac{\pi i n^2}{N}}$,$B_n=e^{\frac{\pi in^2}{N}}$,则原来的傅里叶变换公式可以看作:

$$ X_k=B_n^{*}\sum_{n=0}^{N-1}A_nB_{n-k} $$

后面的式子可以看作是两个序列$A_n,B_n$的卷积结果,因此也就可以将DFT计算转换成卷积计算。

卷积定理

卷积是两个函数或者序列的数学运算,广泛运用于信号处理以及图像处理等领域。卷积定理是傅里叶分析的核心概念之一, 它表明两个信号在时间域的卷积等价于它们在频域的乘积。

假设有两个离散序列$a_n$和$b_n$,用$\star$表示两个序列的卷积,用$F$表示傅里叶变换,用$F^{-1}$表示傅里叶逆变换,则卷积定理可以表示为:

$$ (a_n \star b_n)=F^{-1}(F(a_n) \cdot F(b_n)) $$

利用卷积定理可以将时间域中的离散信号的卷积转换为频域中的点积操作,利用快速傅里叶变换,将计算复杂度从$O(N^2)$降低到$O(NlogN)$。

零填充

FFT计算假设输入信号是周期性的,在使用FFT计算卷积时,通常需要通过零填充来扩展卷积序列以避免循环卷积的影响(因为我们计算的通常是线性卷积),来保证结果的正确性。在普通的卷积计算中,通常会将序列扩展为$N_1+N_2-1$,其中$N_1$和$N_2$分别表示输入序列和卷积序列的长度。对于bluestein算法中的FFT序列来说,假设该FFT序列长度为$N$,则应当扩展为$2N-1$。又因为长度为$2^n$的FFT序列在计算时更为高效,因此在bluestein算法中通常将序列扩充为离$2N-1$最近的2的幂。

实现流程

根据上面算法的原理介绍,具体的实现过程就十分清楚了:

  1. 首先对输入序列$A_n = x_ne^{-\frac{\pi j}{N}n^2}$进行零填充;
  2. 接着对输入序列$B_n=e^{\frac{\pi j}{N}n^2}$进行零填充;
  3. 对执行零填充后的序列分别执行快速傅里叶变换;
  4. 对执行傅里叶变换后的结果执行乘积操作;