傅里叶变换与离散傅里叶变换DFT

一个时域上的函数 $x(t)$ 其傅里叶变换定义为

\tilde{X}(\omega)=\int_{-\infty}^{+\infty}x(t)e^{-i\omega t}d t

如果用频率表示, 将 $\omega=2\pi f$ 代入可得

\tilde{X}(f)=\int_{-\infty}^{+\infty}x(t)e^{-i2\pi f t}d t

对于计算机而言, 记录的下来的函数是对真实函数的采样. 如果记采样的时间周期为 $T$ , 那么采样得到的离散函数自变量应当是采样轮次

x[n]=x(nT)

为了得到离散傅里叶变换, 需要具体地描述采样的过程. 利用 $\delta$ 函数

\delta_s(t)=\sum_{-\infty}^{+\infty}\delta(t-nT)

将它与 $x(t)$ 相乘得到采样后的信号

x_s(t)=\sum_{-\infty}^{+\infty}x(t)\delta(t-nT)

它在除了采样点 $nT$ 之外的点都为零, 虽然在采样点取值正无穷, 但是鉴于

\int_{-\infty}^{+\infty}x(t)\delta(nT)=x(nT)

如果认为采样的过程是对一段段的时间积分,即

x[n]=\int_{nT-T/2}^{nT+T/2}x(t)dt

那么确实有

x[n]=\int_{nT-T/2}^{nT+T/2}x_s(t)dt

当 $nT$ 包含在积分区间内时,积分就为 $x(nT)$ ,否则为零. 因而认为 $x_s(t)$ 是 $x[n]$ 在连续下的形式, 对 $x[n]$ 做傅里叶变换就是在对 $x_s(t)$ 做傅里叶变换. 鉴于 $\delta$ 函数性质优良, 那么

X_s(\omega)=\int_{-\infty}^{+\infty}\left(\sum_{-\infty}^{+\infty}x(t)\delta(t-nT)\right)e^{-i\omega t}dt=\sum_{-\infty}^{+\infty}x(nT)e^{-i\omega nT}

不过显然不可能进行无穷多次采样. 假定在 $0\sim t_0$ 的时间内进行 $N$ 次采样,采样间隔为 $T$ , 其余未采样的地方视为周期的, 即将 $x(t)$ 视为周期函数

x(t)=\begin{cases} x(t)&t\in[0,t_0)\\ x(t+kt_0)&\text{其他},k\in \mathbb{N} \end{cases}

那么它可以进行傅里叶级数展开

x(t)=\sum_{-\infty}^{+\infty}e^{i2\pi n t/t_0}\dfrac{1}{t_0}\int_0^{t_0}x(t)e^{-i2\pi nt/t_0}dt

这说明它的傅里叶变换是离散的

X(\omega)=\sum_{-\infty}^{+\infty}2\pi\delta(\omega-n\omega_0)\dfrac{1}{t_0}\int_0^{t_0}x(t)e^{-in\omega_0 t}dt,\quad \omega_0=\dfrac{2\pi}{t_0}

利用 $\omega=2\pi f$ 将其写为频率

X(f)=\sum_{-\infty}^{+\infty}\delta(f-nf)f_0\int_0^{t_0}x(t)e^{-i2\pi nf_0 t}dt,\quad f_0=\dfrac{1}{t_0}

将周期的 $x_s(t)$ 代入后式得到

f_0\int_0^{t_0}x(t)e^{-i2\pi nf_0 t}dt=f_0\sum_{k=0}^{N-1}x(kT)e^{-i2\pi n kTf_0}

鉴于 $f_0=1/t_0=1/NT$ ,那么

f_0\int_0^{t_0}x(t)e^{-i2\pi nf_0 t}dt=f_0\sum_{k=0}^{N-1}x(kT)e^{-i2\pi n k/N}

记

离散傅里叶变换DFT

X[n]=\sum_{k=0}^{N-1}x[k]e^{-i2\pi n k/N}

鉴于 $e^{-2\pi nk/N}$ 的周期性, 只需要计算 $0<n<N-1$ 即可

FFT的数学基础

根据DFT的表达式不难发现, 如果进行朴素的乘法, 求解一个 $X[n]$ 复杂度为 $O(n)$ , 求解全部 $N$ 个的复杂度为 $O(n^2)$

不过稍加分析,不难发现求解的主要复杂度在于后面的相位项

W[n,k]=\exp\left\{-i\dfrac{2\pi}{N}nk\right\}

由于 $k$ 与 $x[k]$ 相联系, 从 $n$ 处想办法. 注意到

W[n+N/2,k]=\exp\left\{-i\dfrac{2\pi}{N}nk\right\}\exp\left\{-i\pi k\right\}

这意味着当 $k$ 为偶数时有

W[n+N/2,k]=W[n,k]

因而计算 $X[n]$ 时考虑将 $n$ 分奇偶

X[n]=\sum_{k=0}^{N-1}x[k]W[n,k]=\sum_{k=0}^{N/2-1}x[2k]W[n,2k]+\sum_{k=0}^{N/2-1}x[2k+1]W[n,2k+1]

后面的 $2k+1$ 不是偶数, 就把它变成偶数, 利用

W[n,2k+1]=W[n,1]W[n,2k]

那么

X[n]=\sum_{k=0}^{N/2-1}x[2k]W[n,2k]+W[n,1]\sum_{k=0}^{N/2-1}x[2k+1]W[n,2k]

鉴于 $W[n+N/2,k]=W[n,k]$ , 那么

X[n+N/2]=\sum_{k=0}^{N/2-1}x[2k]W[n,2k]+W[n+N/2,1]\sum_{k=0}^{N/2-1}x[2k+1]W[n,2k]

不好处理的是 $W[n+N/2,1]$ , 不过不难发现

W[n+N/2,1]=\exp\left\{-i\dfrac{2\pi}{N}n-\pi i\right\}=-\exp\left\{-\dfrac{2\pi}{N}n\right\}=-W[n,1]

这就意味着

递推式1

X[n+N/2]=X[n]\text{的偶数部分}-W[n,1]\cdot X[n]\text{的奇数部分}

这告诉我们, 知道了 $X[n]$ 的偶数部分就等于知道了 $X[n+N/2]$ 的偶数部分, 奇数也一样;

下一步是求解奇数和偶数部分. 先考察偶数部分

\sum_{k=0}^{N/2-1}x[2k]W[n,2k]

记 $N'=N/2$ , 将偶数位置单独挑出来合成一份 $N'$ 长度的序列, 得到 $x'[k],k=0\sim N'$ , 那么

W[n,2k]=\exp\left\{-i\dfrac{2\pi}{N}n2k\right\}=\exp\left\{-i\dfrac{2\pi}{N'}nk\right\}=W'[n,k]

那么就可以将偶数部分改写为

X[n]\text{的偶数部分}=\sum_{k=0}^{N'-1}x'[k]W'[n,k]

这不正好是偶数序列 $x'[k]$ 的傅里叶变换. 再考察奇数部分

\sum_{k=0}^{N/2-1}x[2k+1]W[n,2k+1]

如法炮制,令 $N'=N/2$ , 将奇数位置挑出来合成 $x'[k]$ , 那么

W[n,2k+1]=W[n,1]W[n,2k]=W[n,1]W'[n,k]

那么就可以将奇数部分改写为

X[n]\text{的奇数部分}=W[n,1]\sum_{k=0}^{N'-1}x'[k]W'[n,k]

进而得到奇数部分和偶数部分的递推式

递推式2

X[n]\text{的偶数部分}=DFT\{x的偶数序列\}[n]

X[n]\text{的奇数部分}=W[n,1]DFT\{x的奇数序列\}[n]

快速傅里叶变换FFT算法

初始化长度为 $N$ 的序列 $x^{(0)}$
对于当前序列 $x^{(m)}$ , 求解它的FFT只需要求解前半部分的奇数和偶数部分, 后半部分即可由递推式自然得到 $X[n+N/2]=X[n]\text{的偶数部分}-W[n,1]\cdot X[n]\text{的奇数部分}$
求解前半部分 $n=0\sim N/2-1$ , 利用递推式 $X[n]\text{的偶数部分}=DFT\{x的偶数序列\}[n]$ $X[n]\text{的奇数部分}=W[n,1]DFT\{x的奇数序列\}[n]$ 问题转化为分别求解奇偶序列的FFT

如此, FFT是一种二分法, 其时间复杂度为 $O(n\log n)$

快速逆傅里叶变换IFFT

注意到

x[n]=\dfrac{1}{N}\sum_{k',k=0}^{N-1}x[k]e^{-i2\pi n(k-k')/N}=\dfrac{1}{N}\sum_{k'=0}^{N-1}X[k']e^{i2\pi nk'/N}

因而定义离散傅里叶逆变换IDFT

离散傅里叶逆变换IDFT

x[n]=\dfrac{1}{N}\sum_{k=0}^{N-1}X[k]e^{i2\pi nk/N}

它长得跟DFT一样,只不过将 $i$ 改为了 $-i$ ,并除了个 $N$ , 那么计算IFFT实际上只需要

IFFT\{X\}=\dfrac{1}{N}FFT\{X;i\to -i\}

目录

傅里叶变换与离散傅里叶变换DFT

FFT的数学基础

快速傅里叶变换FFT算法

快速逆傅里叶变换IFFT