Fourier transform笔记

https://www.youtube.com/watch?v=spUNpyF58BY
笔记

Fourier transform

对于一个声音，可以把它拆解为多个正弦波，上面是两个正弦波组成的声波，但是当情况更复杂以后

很难，从声波中提取出有用的信息，所以采用下面的方法

可以看到，将时间上的波形转换为了左下角的环形的波形，其中，与远点的距离与上面图中与时间轴相对应（白色箭头），图中，有两个频率，一个是时间上的频率，一秒内有3个beats，另一个是打包上图得到左下图的频率（图中虚线），这里是0.5c/s，也就是左下图旋转一周，时间上过了2s。不同的频率，得到的图是不同的。


把左下图中的图形当作一个整体，得到质心，绘制出下图，其中右下图中的y轴是左下图质心的x轴坐标信息。


这样就得到了频域的信息，可以看到，在3明显突起，更改幅度，依然如此

3正好是频率的值，这不是巧合，可以用它来分解波形

为什么是almost，是因为需要去掉高频信息


反向fourier transform，与正向的没有区别

现在已经有了x轴上的fourier变换，需要考虑y轴，既然是二维，可以用复数域来解决。

欧拉公式


$$e^{2\pi ift}$$
其中，f是频率，t是时间，逆时针旋转，那么顺时针旋转
$$e^{?2\pi ift}$$
再乘以强度g（t），这就是我们需要的

这样就得到了fourier变换公式，但多了$\frac{1}{t_2?t_1}$,去掉这一部分就可以表示随着时间变化，质心与原点距离是变化的，不是一成不变的。

laplace transform

laplace变换就是对$e^{?t}sint$的fourier变换
$$\begin{gathered} X(\omega )={\int }_{?\infty }^{\infty }x(t)e^{?i\omega t}dt \\ X(\omega )={\int }_0^{\infty }{e}^{?t}sint{e}^{?i\omega t}dt \\ X(\omega )=\frac{1}{1+(1+i\omega {)}^2} \\ X(\omega )=\frac{1}{2?{\omega }^2+2\omega i} \end{gathered}$$
从第二个式子到第三个式子的推导如下
$$\begin{array}{rl}X(\omega )& ={\int }_0^{\infty }{e}^{?t}sint{e}^{?i\omega t}dt\\ & =?\frac{1}{1+i\omega }{\int }_0^{\infty }sint{e}^{?(1+i\omega )t}d(?(1+i\omega )t)\\ & =?\frac{1}{1+i\omega }sint{e}^{?(1+i\omega )t}{|}_0^{\infty }+?\frac{1}{1+i\omega }{\int }_0^{\infty }sint{e}^{?(1+i\omega )t}dt\\ & =?\frac{1}{1+i\omega }{\int }_0^{\infty }sint{e}^{?(1+i\omega )t}dt\\ & =?\frac{1}{(1+i\omega {)}^2}cost{e}^{?(1+i\omega )t}{|}_0^{\infty }??\frac{1}{(1+i\omega {)}^2}{\int }_0^{\infty }sint{e}^{?(1+i\omega )t}dt\\ & =\frac{1}{(1+i\omega {)}^2}??\frac{1}{(1+i\omega {)}^2}{\int }_0^{\infty }sint{e}^{?(1+i\omega )t}dt\end{array}\text{\hspace{1em}?}$$

很明显，$\omega$控制正弦频率，$\alpha$控制指数衰减

证明如下：
$$\begin{array}{rl}{\int }_0^{\infty }{x}^{\prime }(t)e^{?st}dt& =x(t)e^{?st}{|}_0^{\infty }+s{\int }_0^{\infty }{e}^{?st}x(t)dt\\ & =sx(s)?x(0)\end{array}$$
所以laplace变换可以用来求解常微分方程（ODEs）

FFT(快速傅里叶变换)

FFT是DFT的计算方法，首先了解DFT

DFT（离散傅里叶变换）

实际上就是将积分表示为离散的加和形式
$$X[k]=\sum _{n=0}^{N?1}{x}_n?e^{\frac{?i2\pi }{N}kn}$$
其中，k是频率，总共有N个采样节点，n代表第n个点。
表示为矩阵形式就是：
[ X 0 , X 1 , X 2 , ? ? , X N ? 1 ] ? 1 × N = [ x 0 , x 1 , x 2 , ? ? , x N ? 1 ] ? 1 × N [ k = 0 k = 1 k = 2 ? k = N ? 1 ↓ ↓ ↓ ↓ ] ? N × N \overbrace{[X_0, X_1, X_2, \cdots, X_{N-1}]}^{1\times N} = \overbrace{[x_0, x_1, x_2, \cdots, x_{N-1}]}^{1\times N} \overbrace{\begin{bmatrix} k=0 & k=1 & k=2 & \cdots & k=N-1\\ \Bigg\downarrow & \Bigg\downarrow & \Bigg\downarrow & & \Bigg\downarrow \end{bmatrix}}^{N \times N} [X0?,X1?,X2?,?,XN?1?] ?1×N?=[x0?,x1?,x2?,?,xN?1?] ?1×N? ?k=0↓ ???k=1↓ ???k=2↓ ?????k=N?1↓ ??? ? ?N×N?
如果N很小的时候，这个矩阵乘法是容易计算的，但是现实中，对于信号常常有成百上千的采样点，那么计算量就非常大了，所以就需要FFT。
举例：有8个采样点，用不同频率的正弦函数去得到对应频率的值，每一个值是一个bin

下面需要理解DFT为什么这么做，以及它面临的一些问题

采样频率

首先，根据nyquist定理，要采样频率为f的信号，需要2f+1的采样点，人耳能听到的最高频声音是20khz，所以现实中采样率一般为44.1khz

理想情况

假设有一个f=2hz的余弦，那么我就希望在频域中2的位置有峰值，其他位置为0

那就将离散化得到的向量与不同频率下的波形去求相似性，如果相似度高，就代表该频率下有峰值，这也就是为什么DFT是有效的，意味着采样率一定等于最大频率（N=N）。
但是，通过计算不同频率的波的相似性，会发现问题，$\overrightarrow{a_1}=\overrightarrow{a_7}$,$\overrightarrow{a_2}=\overrightarrow{a_6}$等，导致右下图中很多值不为0

也就是A矩阵不仅不正交，还不可逆，所以不是DFT的最终版本
但是可以发现右下图的左上角这一部分是可逆的，所以每次都采样双倍，但只是采用其中的前半部分。

解决相的问题

通过之前的解决方案，无论是改变幅度，还是频率，或者是初始值，都可以实现FT的功能，但是，当改变相位的时候，频域中的值受到影响，这是不正确的。这是因为在这个变化过程中，无论是sin还是cos都会存在为0的情况，所以要将两者组合在一起。

真正的DFT

就是为sin和cos分别保存一个矩阵，合并后，就得到了我们开始的时候DFT的公式。

FFT引入原因

利用避免重复计算来提升计算速度，上图需要进行64个复数域乘法，改进后，只需要计算24个

FFT原理

对DFT中的部分用W表示，方便表达

简单来讲就是数学家利用上面提到的旋转因子W的周期性，对称性等性质进行公式化简。在算法编程中则是不断利用已经计算过的点来算新的点，即：旧点算新点。
FFT需要依赖蝶形图进行计算
$$\begin{array}{rl}X(k)& =\sum _{n=0}^{N?1}x(n)W_N^{nk}\\ & =\sum _{n=0,2,4,?}^{N?1}x(n)W_N^{nk}+\sum _{n=1,3,5?}^{N?1}x(n)W_N^{nk}\stackrel{\text{令n=2r}}{\to }\\ & =\sum _{r=0,1,?}^{\frac{N}{2}?1}x(2r)W_N^{2rk}+\sum _{r=0,1,?}^{\frac{N}{2}?1}x(2r+1)W_N^{(2r+1)k}\stackrel{x(2r)=x_1(r)}{\to }\\ & =\sum _{r=0,1,?}^{\frac{N}{2}?1}{x}_1(r)W_N^{2rk}+\sum _{r=0,1,?}^{\frac{N}{2}?1}{x}_2(r)W_N^{(2r+1)k}\underset{W_N^{(2r+1)k}=W_N^k?W_{N/2}^{rk}}{\overset{W_N^{2rk}=W_{N/2}^{rk}}{\to }}\\ & =\sum _{r=0,1,?}^{\frac{N}{2}?1}{x}_1(r)W_{N/2}^{rk}+W_N^k\sum _{r=0,1,?}^{\frac{N}{2}?1}{x}_2(r)W_{N/2}^{(rk)}\\ & =X_1(k)+W_N^k{X}_2(k)\end{array}$$
上面的过程没有减少运算量继续进行简化
$$\begin{gathered} \because X_{(}k)=\sum _{r=0,1,?}^{\frac{N}{2}?1}{x}_1(r)W_{N/2}^{rk} \\ \therefore X_1(k+N/2)=\sum _{r=0,1,?}^{\frac{N}{2}?1}{x}_1(r)W_{N/2}^{r(k+N/2)}=\sum _{r=0,1,?}^{\frac{N}{2}?1}{x}_1(r)W_{N/2}^{rk}=X_1(k) \end{gathered}$$
同理，$X_2(k+N/2)=X_2(k)$
$$\begin{gathered} \because X(k)=X_1(k)+W_N^k{X}_2(k) \\ \therefore \begin{array}{rl}X(k+N/2)& =X_1(k+N/2)+W_N^{k+N/2}{X}_2(k+N/2)\stackrel{W_N^{k+N/2}=?W_N^k}{\to }\\ & =X_1(k+N/2)+W_N^k{X}_2(k+N/2)=X_1(k)?W_N^k{X}_2(k)\end{array} \end{gathered}$$
综上推导得到：
$$\begin{gathered} X(k)=X_1(k)+W_N^k{X}_2(k) \\ X(k+N/2)=X_1(k)?W_N^k{X}_2(k) \end{gathered}$$
这里都是用X1和X2来计算，所以缩减了计算量
继续推导


最后一步