离散傅里叶变换

前言

由于实际中接触的信号是离散的，研究离散傅里叶变换是非常必要的。与连续傅里叶变换不同，这时我们只能预测信号，没有理论上的保证。一般的采样方法自然是等间隔采样 \[ x_j=f(jT) \] 即使这样，我们也只能选取其中有限信息周期性延拓之后处理。

周期系列的DFT

对于\(y=\{y_j\},y_j=y_{j+n}\)

DFT (Discrete Fourier Transform) 即离散傅里叶变换定义为 \[ F_n(y)=\bar{y}=\{\bar{y_j}\} \] 其中 \[ y_j=\sum_{k=0}^{n-1}y_k\bar{\omega}^{jk},\omega=e^{\frac{2\pi i}{n}} \] 逆变换IDFT为 \[ y_j=\frac{1}{n}\sum_{k=0}^{n-1}\bar{y_k}\omega^{jk} \] 性质都与连续傅里叶变换类似。

快速傅里叶变换(FFT)

事实上，DFT在很长一段时间内并没有受到人们的重视，直到FFT算法的出现使得DFT的运算大大简化，这才在实际中得到广泛运用。

DFT：\(O(n^2)\) FFT：\(O(n\log n)\)

思想：将序列\(x(j)\)分为两组 \[ \begin{cases} x(2r)=x_1(r)\\ x(2r+1)=x_2(r) \end{cases} \]

\[ \begin{aligned} \bar{x}(k)=\sum_{j=0}^{n-1}x(j)\omega_n^{-jk}&=\sum_{j=0}^{n/2-1}x(2r)\omega_n^{-2rk}+\sum_{j=0}^{n/2-1}x(2r+1)\omega_n^{-(2r+1)k}\\ &=\sum_{j=0}^{n/2-1}x(2r)\omega_n^{-2rk}+\omega_n^{-k}\sum_{j=0}^{n/2-1}x(2r+1)\omega_n^{-2rk}\\ &=G(k)+\omega_n^{-k} H(k) \end{aligned} \]

同时 \[ \bar{x}(k+\frac{n}{2})=G(k)-\omega_n^{-k}H(k) \] 可以看出一个规模为n的离散傅里叶变换问题被划为了两个规模为n/2的子问题，合并需\(O(n)\)时间，故总时间复杂度\(O(n\log n)\)。

如图显示了一次FFT大致的过程。

如果规模不是n的幂可以考虑补零操作

自然顺序	二进制表示	码位倒置	码位倒置顺序
0	000	000	0
1	001	100	4
2	010	010	2
3	011	110	6
4	100	001	1
5	101	101	5
6	110	011	3
7	111	111	7

实际操作中，一般先按自然顺序输入存储单元，然后再通过变址运算将自然顺序的存储转换成码位倒置顺序的存储，然后进行FFT的原位计算。

一个小应用——切比雪夫插值

Chebyshev多项式

\[ \begin{aligned} T_n(\cos\theta)&=\cos(n\theta)\\ T_0(1&)=1\\ T_1(1&)=x\\ T_2(1)=&2x^2-1\\ T_3(1)=&4x^3-3x\\ \cdots \end{aligned} \]

插值节点选在\(\cos k\pi/n\)上，这样可以尽量使误差变小，一定程度上规避龙格现象的产生。（数值分析内容就不展开了）

对一个函数计算插值系数 \[ P(x)=\sum_{j=0}^{n}a_jT_j(x) \] 令\(x_k=\cos(k\pi/n)\)

即 \[ y_k=P(x_k)=\sum_{j=0}^{n}a_jT_j(x_k)=\sum_{j=0}^{n}a_j\cos(jk\pi/n) \] 则可以改写为 \[ y_k=\frac{1}{2}\sum_{j=0}^na_j\omega^{jk}+\frac{1}{2}\sum_{j=-n}^{0}a_{-j} \omega^{jk}=\sum_{j=-n}^nc_j\omega^{jk} \] 将\((y_k)\)延拓成\(2n\)周期的偶序列则有 \[ y_k=\sum_{j=-n}^nc_j\omega^{jk},0\le j\le 2n-1 \] 用DFT公式套上去即可得到 \[ a_j=\frac{1}{n\epsilon_j}\sum_{k=0}^{2n-1}y_k\bar{\omega}^{jk},j=0,1,\cdots,n \] 其中 \[ \epsilon_0=\epsilon_n=2,\epsilon_1=\cdots=\epsilon_{n-1}=1 \] 即可利用快速离散傅里叶变换得到系数。

离散滤波器

假设\(X,Y\)为离散信号空间，称算子\(F:X\rightarrow Y\)是线性时不变的，如果满足 \[ \begin{aligned} &F(\alpha x+\beta y)=\alpha F(x)+\beta F(y)\\ &F(T_p(x))=T_p(F(x)) \end{aligned} \] 其中 \[ (T_p(x))_k=x_{k-p} \] 。离散傅里叶变换与连续傅里叶变换也有非常类似的结论：

$ Thm $ 如果\(F\)是离散信号空间上的线性时不变算子，则存在序列\(f\)，使得 \[ F(x)=f*x \]

Z变换

序列\(x\)的Z变换定义为函数\(\hat{x}:[-\pi,\pi]\rightarrow \mathbb{C}\)： \[ \hat{x}(\phi)=\sum_{j\in \mathbb{Z}}x_je^{-ij\phi} \] 注：\(z=e^{i\phi}\)有 \[ \hat{x}(\phi)=\sum_{j\in \mathbb{Z}}x_jz^{-j} \]

Z变换与Fourier级数

假设\(f\in L^2[-\pi,\pi]\)，其Fourier级数展开为 \[ f(\phi)=\sum_{n\in \mathbb{Z}}x_ne^{in\phi} \] 即傅里叶级数是一种\(L^2[-\pi,\pi]\rightarrow l^2\)的变换。

Z变换把序列转化为函数\(f(-\cdot)\in L^2[-\pi,\pi]\)

卷积算子与Z变换

\(Thm\) 假设\(f=(f_n),x=(x_n)\in l^2\)，则 \[ (\widehat{f*x})(\phi)=\hat{f}(\phi)\hat{x}(\phi) \] 验证是容易的。

窗口傅里叶变换（WFT）

使用傅里叶变换确实是信号处理一个非常常用的手法，但是它只能适合平稳信号的分析，给出一个信号总体上包含哪些频率的成分，并不能给出各成分出现的时刻。但是实际上自然中大部分信号是非平稳的，任何一个频谱都是和所有的时间相关。为了更好地分析信号，我们需要考虑将信号分成若干份，每一段分别做傅里叶变换，这就引出了窗口傅里叶变换。

窗函数

非平凡函数\(\omega\in L^2(\mathbb{R})\)称为窗函数如果满足\(t\omega\in L^2(\mathbb{R})\)。

中心 \[ t^*=\frac{1}{\Vert \omega\Vert_{L^2}^2}\int_{\mathbb{R}}t|\omega(t)|^2dt \] 半径 \[ \Delta_\omega=\frac{1}{\Vert \omega\Vert_{L^2}}(\int_{\mathbb{R}}(t-t^*)^2|\omega(t)|^2dt)^{1/2} \]

时间窗、频率窗

设\(g\)是实的窗函数，且\(g(t)=g(-t),\Vert g\Vert=1\)则\(f\)的窗口傅里叶变换定义为 \[ S[f](\lambda,b)=\frac{1}{\sqrt{2\pi}}\int_{\mathbb{R}}f(t)\overline{g(t-b)}e^{-i\lambda t}dt \]

一般\(g\)是一个在某区间上取值，区间之外快速趋于零的函数，例如Guass型函数。表达式中\(f,g\)相乘一项就是选取\(f\)在\(b\)周围一块，再做傅里叶变换。随着\(b\)变化，选取区间（窗口）也在变化。

令\(W_{\lambda,b}(t)=g(t-b)e^{i\lambda t}\)

则由窗口傅里叶变换定义可得 \[ S[f](\lambda,b)=\frac{1}{\sqrt{2\pi}}\int_{\mathbb{R}}f(t)\overline{W_{\lambda,b}(t)}dt \] \(W_{\lambda, b}\)的中心和半径分别是\(t^*+b\)和\(\Delta_g\)，则窗口傅里叶变换主要给出了\(f(t)\)在时间窗 \[ [t^*+b-\Delta_g,t^*+b+\Delta_g] \] 里面的局部信息。

另一方面，令 \[ V_{\lambda,b}(\omega)=\mathfrak{F}[W_{\lambda,b}](\omega)=e^{i\lambda b}e^{-i\omega b}\hat{g}(\omega-\lambda) \] \(V_{\lambda, b}\)的中心和半径分别是\(\lambda^*+\lambda\)和\(\Delta_\hat{g}\) \[ S[f](\lambda,b)=\frac{1}{\sqrt{2\pi}}\langle f,W_{\lambda,b}(t)\rangle=\frac{1}{\sqrt{2\pi}}\langle \hat{f},V_{\lambda,b}(t)\rangle \] 即窗口傅里叶变换还给出了频率窗 \[ [\lambda^*+\lambda-\Delta_{\hat{g}},\lambda^*+\lambda+\Delta_{\hat{g}}] \] 里面的局部化信息。

可以观察到时频指标移动时，窗函数的时频窗覆盖了整个时频平面。可以证明，一个信号可以用它的窗口傅里叶变换来恢复。

但是时频窗的总面积是固定的，对于窗口宽度的确定就成了一个非常令人头疼的课题。窗太窄，窗内的信号太短，会导致频率分析不够精准，频率分辨率差。窗太宽，时域上又不够精细，时间分辨率低。对于该问题的解决方案就引出了小波变换。

（傅里叶变换到此结束，进入分类主题——小波分析）