Fourier transform笔记

2023-12-21 18:37:44

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

Fourier transform

图1
对于一个声音,可以把它拆解为多个正弦波,上面是两个正弦波组成的声波,但是当情况更复杂以后
在这里插入图片描述
很难,从声波中提取出有用的信息,所以采用下面的方法
在这里插入图片描述
可以看到,将时间上的波形转换为了左下角的环形的波形,其中,与远点的距离与上面图中与时间轴相对应(白色箭头),图中,有两个频率,一个是时间上的频率,一秒内有3个beats,另一个是打包上图得到左下图的频率(图中虚线),这里是0.5c/s,也就是左下图旋转一周,时间上过了2s。不同的频率,得到的图是不同的。
在这里插入图片描述
在这里插入图片描述
把左下图中的图形当作一个整体,得到质心,绘制出下图,其中右下图中的y轴是左下图质心的x轴坐标信息。
在这里插入图片描述
在这里插入图片描述
这样就得到了频域的信息,可以看到,在3明显突起,更改幅度,依然如此
在这里插入图片描述
3正好是频率的值,这不是巧合,可以用它来分解波形
在这里插入图片描述
为什么是almost,是因为需要去掉高频信息
在这里插入图片描述
在这里插入图片描述
反向fourier transform,与正向的没有区别
在这里插入图片描述
现在已经有了x轴上的fourier变换,需要考虑y轴,既然是二维,可以用复数域来解决。
在这里插入图片描述
欧拉公式
在这里插入图片描述
在这里插入图片描述
e 2 π i f t e^{2\pi ift} e2πift
其中,f是频率,t是时间,逆时针旋转,那么顺时针旋转
e ? 2 π i f t e^{-2\pi ift} e?2πift
再乘以强度g(t),这就是我们需要的
在这里插入图片描述
这样就得到了fourier变换公式,但多了 1 t 2 ? t 1 \frac{1}{t_2-t_1} t2??t1?1?,去掉这一部分就可以表示随着时间变化,质心与原点距离是变化的,不是一成不变的。
在这里插入图片描述

laplace transform

在这里插入图片描述
laplace变换就是对 e ? t s i n t e^{-t}sint e?tsint的fourier变换
X ( ω ) = ∫ ? ∞ ∞ x ( t ) e ? i ω t d t X ( ω ) = ∫ 0 ∞ e ? t s i n t e ? i ω t d t X ( ω ) = 1 1 + ( 1 + i ω ) 2 X ( ω ) = 1 2 ? ω 2 + 2 ω i X(\omega) = \int^{\infty}_{-\infty}x(t)e^{-i\omega t}dt \\ X(\omega) = \int^{\infty}_{0}e^{-t}sinte^{-i\omega t} dt \\ X(\omega) = \frac{1}{1+(1+i\omega)^2} \\ X(\omega) = \frac{1}{2-\omega^2 + 2\omega i} X(ω)=??x(t)e?tdtX(ω)=0?e?tsinte?tdtX(ω)=1+(1+)21?X(ω)=2?ω2+2ωi1?
从第二个式子到第三个式子的推导如下
X ( ω ) = ∫ 0 ∞ e ? t s i n t e ? i ω t d t = ? 1 1 + i ω ∫ 0 ∞ s i n t e ? ( 1 + i ω ) t d ( ? ( 1 + i ω ) t ) = ? 1 1 + i ω s i n t e ? ( 1 + i ω ) t ∣ 0 ∞ + ? 1 1 + i ω ∫ 0 ∞ s i n t e ? ( 1 + i ω ) t d t = ? 1 1 + i ω ∫ 0 ∞ s i n t e ? ( 1 + i ω ) t d t = ? 1 ( 1 + i ω ) 2 c o s t e ? ( 1 + i ω ) t ∣ 0 ∞ ? ? 1 ( 1 + i ω ) 2 ∫ 0 ∞ s i n t e ? ( 1 + i ω ) t d t = 1 ( 1 + i ω ) 2 ? ? 1 ( 1 + i ω ) 2 ∫ 0 ∞ s i n t e ? ( 1 + i ω ) t d t \begin{align} X(\omega) &= \int^{\infty}_{0}e^{-t}sinte^{-i\omega t}dt \nonumber \\ &= -\frac{1}{1+i\omega}\int^{\infty}_{0}sinte^{-(1+i\omega) t}d(-(1+i\omega) t) \nonumber \\ &= -\frac{1}{1+i\omega}sint e^{-(1+i\omega) t} |^{\infty}_{0}+ -\frac{1}{1+i\omega}\int^{\infty}_{0}sinte^{-(1+i\omega) t}dt \nonumber \\ &= -\frac{1}{1+i\omega}\int^{\infty}_{0}sinte^{-(1+i\omega) t}dt \nonumber \\ &= -\frac{1}{(1+i\omega)^2}cost e^{-(1+i\omega) t} |^{\infty}_{0} - -\frac{1}{(1+i\omega)^2}\int^{\infty}_{0}sinte^{-(1+i\omega) t}dt\nonumber \\ &= \frac{1}{(1+i\omega)^2}- -\frac{1}{(1+i\omega)^2}\int^{\infty}_{0}sinte^{-(1+i\omega) t}dt\nonumber \\ \end{align} X(ω)?=0?e?tsinte?tdt=?1+1?0?sinte?(1+)td(?(1+)t)=?1+1?sinte?(1+)t0?+?1+1?0?sinte?(1+)tdt=?1+1?0?sinte?(1+)tdt=?(1+)21?coste?(1+)t0???(1+)21?0?sinte?(1+)tdt=(1+)21???(1+)21?0?sinte?(1+)tdt??
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

在这里插入图片描述
很明显, ω \omega ω控制正弦频率, α \alpha α控制指数衰减
在这里插入图片描述
证明如下:
∫ 0 ∞ x ′ ( t ) e ? s t d t = x ( t ) e ? s t ∣ 0 ∞ + s ∫ 0 ∞ e ? s t x ( t ) d t = s x ( s ) ? x ( 0 ) \begin{align} \int^{\infty}_{0}x'(t)e^{-st}dt&=x(t)e^{-st}|^{\infty}_{0} + s\int^{\infty}_{0}e^{-st}x(t)dt \nonumber \\ &= sx(s)-x(0) \nonumber \end{align} 0?x(t)e?stdt?=x(t)e?st0?+s0?e?stx(t)dt=sx(s)?x(0)?
所以laplace变换可以用来求解常微分方程(ODEs)

FFT(快速傅里叶变换)

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

DFT(离散傅里叶变换)

实际上就是将积分表示为离散的加和形式
X [ k ] = ∑ n = 0 N ? 1 x n ? e ? i 2 π N k n X[k] = \sum^{N-1}_{n=0}x_n\cdot e^{\frac{-i2\pi}{N}kn} X[k]=n=0N?1?xn??eN?i2π?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)。
但是,通过计算不同频率的波的相似性,会发现问题, a 1 ? = a 7 ? \vec{a_1}=\vec{a_7} a1? ?=a7? ?, a 2 ? = a 6 ? \vec{a_2}=\vec{a_6} a2? ?=a6? ?等,导致右下图中很多值不为0
在这里插入图片描述
也就是A矩阵不仅不正交,还不可逆,所以不是DFT的最终版本
但是可以发现右下图的左上角这一部分是可逆的,所以每次都采样双倍,但只是采用其中的前半部分。

解决相的问题

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

真正的DFT

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

FFT引入原因

利用避免重复计算来提升计算速度,上图需要进行64个复数域乘法,改进后,只需要计算24个
在这里插入图片描述

FFT原理

对DFT中的部分用W表示,方便表达
在这里插入图片描述
简单来讲就是数学家利用上面提到的旋转因子W的周期性,对称性等性质进行公式化简。在算法编程中则是不断利用已经计算过的点来算新的点,即:旧点算新点。
FFT需要依赖蝶形图进行计算
X ( k ) = ∑ n = 0 N ? 1 x ( n ) W N n k = ∑ n = 0 , 2 , 4 , ? N ? 1 x ( n ) W N n k + ∑ n = 1 , 3 , 5 ? N ? 1 x ( n ) W N n k → 令n=2r = ∑ r = 0 , 1 , ? N 2 ? 1 x ( 2 r ) W N 2 r k + ∑ r = 0 , 1 , ? N 2 ? 1 x ( 2 r + 1 ) W N ( 2 r + 1 ) k → x ( 2 r ) = x 1 ( r ) = ∑ r = 0 , 1 , ? N 2 ? 1 x 1 ( r ) W N 2 r k + ∑ r = 0 , 1 , ? N 2 ? 1 x 2 ( r ) W N ( 2 r + 1 ) k → W N ( 2 r + 1 ) k = W N k ? W N / 2 r k W N 2 r k = W N / 2 r k = ∑ r = 0 , 1 , ? N 2 ? 1 x 1 ( r ) W N / 2 r k + W N k ∑ r = 0 , 1 , ? N 2 ? 1 x 2 ( r ) W N / 2 ( r k ) = X 1 ( k ) + W N k X 2 ( k ) \begin{align} X(k) &= \sum^{N-1}_{n=0}x(n)W_N^{nk} \nonumber \\ &= \sum^{N-1}_{n=0,2,4,\cdots}x(n)W_N^{nk} + \sum^{N-1}_{n=1, 3, 5\cdots}x(n)W_N^{nk} \xrightarrow{\text{令n=2r}}\nonumber \\ &= \sum^{\frac{N}{2}-1}_{r=0,1,\cdots}x(2r)W_N^{2rk} + \sum^{\frac{N}{2}-1}_{r=0,1, \cdots}x(2r+1)W_N^{(2r+1)k} \xrightarrow{x(2r)=x_1(r)}\nonumber \\ &= \sum^{\frac{N}{2}-1}_{r=0,1,\cdots}x_1(r)W_N^{2rk} + \sum^{\frac{N}{2}-1}_{r=0,1, \cdots}x_2(r)W_N^{(2r+1)k} \xrightarrow[W_N^{(2r+1)k} = W_{N}^{k} * W_{N/2}^{rk}]{W_N^{2rk}=W_{N/2}^{rk}}\nonumber \\ &= \sum^{\frac{N}{2}-1}_{r=0,1,\cdots}x_1(r)W_{N/2}^{rk} + W_N^k\sum^{\frac{N}{2}-1}_{r=0,1, \cdots}x_2(r)W_{N/2}^{(rk)} \nonumber \\ &=X_1(k)+W_N^kX_2(k) \nonumber \end{align} X(k)?=n=0N?1?x(n)WNnk?=n=0,2,4,?N?1?x(n)WNnk?+n=1,3,5?N?1?x(n)WNnk?n=2r ?=r=0,1,?2N??1?x(2r)WN2rk?+r=0,1,?2N??1?x(2r+1)WN(2r+1)k?x(2r)=x1?(r) ?=r=0,1,?2N??1?x1?(r)WN2rk?+r=0,1,?2N??1?x2?(r)WN(2r+1)k?WN2rk?=WN/2rk? WN(2r+1)k?=WNk??WN/2rk??=r=0,1,?2N??1?x1?(r)WN/2rk?+WNk?r=0,1,?2N??1?x2?(r)WN/2(rk)?=X1?(k)+WNk?X2?(k)?
上面的过程没有减少运算量继续进行简化
∵ X ( k ) = ∑ r = 0 , 1 , ? N 2 ? 1 x 1 ( r ) W N / 2 r k ∴ X 1 ( k + N / 2 ) = ∑ r = 0 , 1 , ? N 2 ? 1 x 1 ( r ) W N / 2 r ( k + N / 2 ) = ∑ r = 0 , 1 , ? N 2 ? 1 x 1 ( r ) W N / 2 r k = X 1 ( k ) \because X_(k) = \sum^{\frac{N}{2}-1}_{r=0,1,\cdots}x_1(r)W_{N/2}^{rk} \\ \therefore X_1(k+N/2) = \sum^{\frac{N}{2}-1}_{r=0,1,\cdots}x_1(r)W_{N/2}^{r(k+N/2)}=\sum^{\frac{N}{2}-1}_{r=0,1,\cdots}x_1(r)W_{N/2}^{rk} = X_1(k) X(?k)=r=0,1,?2N??1?x1?(r)WN/2rk?X1?(k+N/2)=r=0,1,?2N??1?x1?(r)WN/2r(k+N/2)?=r=0,1,?2N??1?x1?(r)WN/2rk?=X1?(k)
同理, X 2 ( k + N / 2 ) = X 2 ( k ) X_2(k+N/2) = X_2(k) X2?(k+N/2)=X2?(k)
∵ X ( k ) = X 1 ( k ) + W N k X 2 ( k ) ∴ X ( k + N / 2 ) = X 1 ( k + N / 2 ) + W N k + N / 2 X 2 ( k + N / 2 ) → W N k + N / 2 = ? W N k = X 1 ( k + N / 2 ) + W N k X 2 ( k + N / 2 ) = X 1 ( k ) ? W N k X 2 ( k ) \because X(k)=X_1(k)+W_N^kX_2(k) \\ \therefore \begin{align} X(k+N/2)&=X_1(k+N/2)+W_N^{k+N/2}X_2(k+N/2) \xrightarrow{W_N^{k+N/2}=-W_N^k} \nonumber \\ &= X_1(k+N/2)+W_N^{k}X_2(k+N/2) = X_1(k)-W_N^kX_2(k) \nonumber \end{align} X(k)=X1?(k)+WNk?X2?(k)X(k+N/2)?=X1?(k+N/2)+WNk+N/2?X2?(k+N/2)WNk+N/2?=?WNk? ?=X1?(k+N/2)+WNk?X2?(k+N/2)=X1?(k)?WNk?X2?(k)?
综上推导得到:
X ( k ) = X 1 ( k ) + W N k X 2 ( k ) X ( k + N / 2 ) = X 1 ( k ) ? W N k X 2 ( k ) X(k)=X_1(k)+W_N^kX_2(k) \\ X(k+N/2) = X_1(k)-W_N^kX_2(k) X(k)=X1?(k)+WNk?X2?(k)X(k+N/2)=X1?(k)?WNk?X2?(k)
这里都是用X1和X2来计算,所以缩减了计算量
继续推导
在这里插入图片描述
在这里插入图片描述
最后一步
在这里插入图片描述
在这里插入图片描述

文章来源:https://blog.csdn.net/qq_42725437/article/details/135121474
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。