coverPiccoverPic

理解傅里叶变换

🎈前言

傅立叶变换(Fourier Transform,FT)可以将时域上的信号转变为频域上的信号,使用者从而获得感兴趣的变量,在音视频、图像、神经影像等信号处理领域均有广泛的应用。而快速傅立叶变换(Fast Fourier Transform,FFT)是一种离散信号的傅立叶变换的的快速计算方法。

傅立叶变换在我之前学习中接触过,基本上每天都在和它打交道,但是一直对它的原理不甚了解,因此趁着工作休息之余,用两篇博客的时间,稍微回顾一下。

在前端开发中,其实也可以看到它的身影,只不过对它的使用就只局限于浏览器提供的AudioContext等少数几个 API 就是了…

傅里叶变换的公式是

F(ω)=f(t)eiωtdtF(\omega) = \int_{-\infin}^{\infin}f(t)e^{-i\omega t}dt

其中 f(t)f(t) 是信号时域上的原函数,tt 是时间,ω\omega 是频率。直观上,就像把原来的函数 f(t)f(t) 从时间到振幅的映射“换元”变成了从频率映射,到组成原信号各个正弦波振幅和相位的函数。

✨以 2π 为周期的傅立叶级数

傅立叶变换基于傅立叶级数的推广。从傅立叶级数(Fourier Series)说起,法国数学家傅立叶认为任何周期函数都可以用正弦函数和余弦函数构成的无穷级数来表示,这里假设 f(t)f(t) 是周期为 2π2\pi 的函数,即:

f(t)=A0+n=1Ansin(nt+φn)=A0+n=1An(cosφnsin(nt)+sinφncos(nt))=A0+Ansinφnn=1cos(nt)+Ancosφnn=1sin(nt)f(t)=A_0+\sum_{n=1}^{\infin}A_n\sin(nt+\varphi_n)\\ =A_0+\sum_{n=1}^{\infin}A_n(\cos\varphi_n\sin(nt)+\sin\varphi_n\cos(nt))\\ =A_0+A_n\sin\varphi_n\sum_{n=1}^{\infin}\cos(nt)+A_n\cos\varphi_n\sum_{n=1}^{\infin}\sin(nt)\\

a0=2Ana_0 = 2A_nan=Ansinφna_n = A_n\sin\varphi_nbn=Ancosφnb_n = A_n\cos\varphi_n

f(t)=a02+n=1(ancos(nt)+bnsin(nt))f(t)=\frac{a_0}{2} + \sum_{n=1}^{\infin}(a_ncos(nt)+ b_nsin(nt))

记为①,把 ana_nbnb_n 求解出了就获得组成原函数的各个正弦波的分量了。要求 a0a_0

ππf(t)dt=a02ππdt+ππn=1(ancos(nt)+bnsin(nt))dt\int_{-\pi}^{\pi}f(t)dt=\frac{a_0}{2}\int_{-\pi}^{\pi}dt + \int_{-\pi}^{\pi}\sum_{n=1}^{\infin}(a_n\cos(nt)+ b_n\sin(nt))dt

n=1(ancos(nt)+bnsin(nt))\sum_{n=1}^{\infin}(a_ncos(nt)+ b_nsin(nt)) 一致收敛的条件下,积分和求和可交换顺序:

ππf(t)dt=a02ππdt+n=1(anππcos(nt)dt+bnππsin(nt)dt)\int_{-\pi}^{\pi}f(t)dt=\frac{a_0}{2}\int_{-\pi}^{\pi}dt + \sum_{n=1}^{\infin}(a_n\int_{-\pi}^{\pi}\cos(nt)dt+ b_n\int_{-\pi}^{\pi}\sin(nt)dt)

ππcos(nt)dt\int_{-\pi}^{\pi}\cos(nt)dtππsin(nt)dt\int_{-\pi}^{\pi}\sin(nt)dtn>0n > 0 时为 0,

a0=1πππf(t)dta_0 = \frac{1}{\pi}\int_{-\pi}^{\pi}f(t)dt

要求 aka_k

f(t)cos(kt)=a02cos(kt)+cos(kt)n=1(ancos(nt)+bnsin(nt))f(t)\cos(kt)=\frac{a_0}{2}\cos(kt) + \cos(kt)\sum_{n=1}^{\infin}(a_n\cos(nt)+ b_n\sin(nt))

ππf(t)cos(kt)dt=ππa02cos(kt)dt+ann=1ππcos(nt)cos(kt)dt+bnn=1ππsin(nt)cos(kt)dt=ann=1ππcos(nt)cos(kt)dt+bnn=1ππsin(nt)cos(kt)dt\int_{-\pi}^{\pi}f(t)\cos(kt)dt=\int_{-\pi}^{\pi}\frac{a_0}{2}\cos(kt)dt + a_n\sum_{n=1}^{\infin}\int_{-\pi}^{\pi}cos(nt)\cos(kt)dt+b_n\sum_{n=1}^{\infin}\int_{-\pi}^{\pi}sin(nt)\cos(kt)dt\\ =a_n\sum_{n=1}^{\infin}\int_{-\pi}^{\pi}cos(nt)\cos(kt)dt+b_n\sum_{n=1}^{\infin}\int_{-\pi}^{\pi}sin(nt)\cos(kt)dt

由三角函数正交性,ππcos(nt)cos(kt)dt\int_{-\pi}^{\pi}cos(nt)\cos(kt)dt,仅在 n=kn=k 时不为零,ππsin(nt)cos(kt)dt=0\int_{-\pi}^{\pi}sin(nt)\cos(kt)dt = 0

ππf(t)cos(kt)dt=akππcos(kt)cos(kt)dt=ak(t/2+14ksin(kt)cos(kt))ππ=πak\int_{-\pi}^{\pi}f(t)\cos(kt)dt=a_k\int_{-\pi}^{\pi}cos(kt)\cos(kt)dt\\ =a_k(t/2 + \frac{1}{4k}\sin(kt)\cos(kt))\lvert^{\pi}_{-\pi}\\ =\pi a_k\\

因此

ak=1πππf(t)cos(kt)dta_k=\frac{1}{\pi}\int_{-\pi}^{\pi}f(t)\cos(kt)dt

可以发现,a0{a_0} 也是包括在 ak{a_k} 这种形式中的,类似地,把①左右两边乘 sin(kt)\sin(kt) 再积分,有:

bk=1πππf(t)sin(kt)dtb_k=\frac{1}{\pi}\int_{-\pi}^{\pi}f(t)\sin(kt)dt

🎀以 2L 为周期的情况

对于周期不一定是 2π2\pi 的周期函数,现在记 f(t)f(t) 的周期是 2L2L,则 g(t)=f(tLπ)g(t) = f(\frac{ tL}{\pi}) 周期为 2π2\pi,有

g(t)=a02+n=1(ancos(nt)+bnsin(ntx))ak=1πππg(t)cos(kt)dtbk=1πππg(t)sin(kt)dtg(t)=\frac{a_0}{2} + \sum_{n=1}^{\infin}(a_n\cos(nt)+ b_n\sin(ntx))\\ a_k=\frac{1}{\pi}\int_{-\pi}^{\pi}g(t)\cos(kt)dt\\ b_k=\frac{1}{\pi}\int_{-\pi}^{\pi}g(t)\sin(kt)dt

t=πxLt=\frac{\pi x}{L},换元可得:

f(x)=a02+n=1(ancos(nπxL)+bnsin(nπxL))ak=1LLLf(x)cos(nπxL)dxbk=1LLLf(x)sin(nπxL)dxf(x)=\frac{a_0}{2} + \sum_{n=1}^{\infin}(a_n\cos(\frac{n\pi x}{L})+ b_n\sin(\frac{n\pi x}{L}))\\ a_k=\frac{1}{L}\int_{-L}^{L}f(x)\cos(\frac{n\pi x}{L})dx\\ b_k=\frac{1}{L}\int_{-L}^{L}f(x)\sin(\frac{n\pi x}{L})dx

🚪复数域上傅立叶展开

由欧拉公式 eix=cosx+isinxe^{ix}=\cos x+i\sin x,有 cosx=12(eix+eix)\cos x=\frac{1}{2}(e^{ix}+e^{-ix})sinx=12i(eixeix)\sin x=-\frac{1}{2}i(e^{ix}-e^{-ix})

f(t)=a02+n=1(ancos(nπtL)+bnsin(nπtL))=a02+n=1(an12(einπtL+einπtL)bn12i(einπtLeinπtL))=a02+n=1(anibn2einπtL)+n=1(an+ibn2einπtL)=n=00a02einπtL+n=1(anibn2einπtL)+n=1(an+ibn2einπtL)f(t)=\frac{a_0}{2} + \sum_{n=1}^{\infin}(a_n\cos(\frac{n\pi t}{L})+ b_n\sin(\frac{n\pi t}{L}))\\ =\frac{a_0}{2} + \sum_{n=1}^{\infin}(a_n\frac{1}{2}(e^{i\frac{n\pi t}{L}}+e^{-i\frac{n\pi t}{L}})- b_n\frac{1}{2}i(e^{i\frac{n\pi t}{L}}-e^{-i\frac{n\pi t}{L}}))\\ =\frac{a_0}{2} + \sum_{n=1}^{\infin}(\frac{a_n-ib_n}{2}e^{i\frac{n\pi t}{L}}) + \sum_{n=1}^{\infin}(\frac{a_n+ib_n}{2}e^{-i\frac{n\pi t}{L}})\\ =\sum_{n=0}^{0}\frac{a_0}{2}e^{i\frac{n\pi t}{L}} + \sum_{n=1}^{\infin}(\frac{a_n-ib_n}{2}e^{i\frac{n\pi t}{L}}) + \sum_{n=-1}^{-\infin}(\frac{a_{-n}+ib_{-n}}{2}e^{i\frac{n\pi t}{L}})\\

可以发现,傅立叶级数可以写成这样子的形式:

f(t)=n=(cneinπt/L)cn=12LLLf(t)einπt/Ldtf(t)=\sum_{n=-\infin}^{\infin}(c_ne^{in\pi t/L})\\ c_n=\frac{1}{2L}\int^{L}_{-L}f(t)e^{-in\pi t/L}dt

ω=π/L\omega = \pi /L

f(t)=n=(cneinωt)cn=12LLLf(t)einωtdtf(t)=\sum_{n=-\infin}^{\infin}(c_ne^{in\omega t})\\ c_n=\frac{1}{2L}\int^{L}_{-L}f(t)e^{-in\omega t}dt

于是,cnc_n 就是周期为 2L2L 的函数各个正弦波分量的系数。

🔍傅立叶变换

实际应用中,信号通常都是一段非周期性的数值,对于非周期的信号,处理方法是把整段信号看成是一个周期,从中,就可以得到傅立叶的表达 F(ω)F(\omega),记 cn=12LLLf(t)einω0tdtc_n=\frac{1}{2L}\int^{L}_{-L}f(t)e^{-in\omega_0 t}dtω=nω0\omega = n\omega_0

F(ω)=limL2Lcn=limLLLf(t)einω0tdt=f(t)eiωtdtF(\omega) = \lim_{L\rightarrow \infin} 2Lc_n =\lim_{L\rightarrow \infin} \int^{L}_{-L}f(t)e^{-in\omega_0 t}dt= \int^{\infin}_{-\infin}f(t)e^{-i\omega t}dt

此时,对 f(t)f(t) 来说,cnc_n可以表示成如下形式:

cn=12LF(nω0)c_n=\frac{1}{2L}F(n\omega_0)

表示原信号的函数可以用 F(ω)F(\omega) 表示,这就是傅立叶变换的逆变换:

f(t)=limLn=(12LF(nπ/L)einπt/L)f(t) =\lim_{L\rightarrow \infin}\sum_{n=-\infin}^{\infin}(\frac{1}{2L}F(n\pi/L)e^{in\pi t/L})\\

根据定积分的定义,LL\rightarrow\infin,上式是一个定积分,

f(t)=12(F(nπ/L)einπt/L)dnLf(t) =\frac{1}{2}\int_{-\infin}^{\infin}(F(n\pi/L)e^{in\pi t/L})d\frac{n}{L}

同样地,记 ω=nπ/L\omega=n\pi/L

f(t)=12πF(ω)eiωtdωf(t)=\frac{1}{2\pi}\int^{\infin}_{-\infin}F(\omega)e^{i\omega t}d\omega

🛠离散数据的傅里叶变换

计算机只能处理离散的数据,所需处理的采样信号可以看作以 TT 为周期的连续的时间信号 f(t)f(t) 与冲激串信号相乘的结果。

在一个周期 [t0,t0+T][t_0,t_0+T] 内以 TsT_s 为间隔,采样 NN 次,冲激串信号定义如下:

k=0K1δ(tkTs)\sum_{k=0}^{K -1}\delta(t -kT_s)

其中 δ\delta 是狄拉克函数,它除了 x=0x=0 处不为 0 外都为 0,并且在整个定义域积分结果为 1,而且它有一个特性:

f(x)δ(xx0)dx=f(x0)\int^{\infin}_{-\infin}f(x)\delta(x-x_0)dx = f(x_0)

总之,这是采样信号:

g(t)=f(t)k=0N1δ(tkTs)g(t) = f(t)\sum_{k=0}^{N -1}\delta(t -kT_s)

傅立叶展开中

cn=1TT/2T/2(g(t)k=0N1δ(tkTs))e2πint/Tdt=1Tk=0N1T/2T/2g(t)δ(tkTs)e2πint/Tdt=1Tk=0N1g(kTs)e2πinkTs/T=1Tk=0N1g(kTs)e2πink/Nc_n=\frac{1}{T}\int^{T/2}_{-T/2}(g(t)\sum_{k=0}^{N -1}\delta(t -kT_s))e^{-2\pi int/T}dt\\ =\frac{1}{T}\sum_{k=0}^{N -1}\int^{T/2}_{-T/2}g(t)\delta(t -kT_s)e^{-2\pi int/T}dt\\ =\frac{1}{T}\sum_{k=0}^{N -1}g(kT_s)e^{-2\pi in kT_s/T}\\ =\frac{1}{T}\sum_{k=0}^{N -1}g(kT_s)e^{-2\pi in k/N}

g(kTs)=f[k]g(kT_s) = f[k]

cn=1Tk=0N1f[k]e2πink/Nc_n=\frac{1}{T}\sum_{k=0}^{N -1}f[k]e^{-2\pi in k/N}

ω=2πn/N\omega = 2\pi n/N,定义周期信号的离散傅里叶变换(Discrete Fourier Transform,DFT)为

F(ω)=Tcn=k=0N1f[k]eiωkF(\omega) =Tc_n=\sum_{k=0}^{N -1}f[k]e^{-i\omega k}

逆变换为

f[k]=1Nk=0N1F(ω)eiωkf[k] =\frac{1}{N}\sum_{k=0}^{N -1}F(\omega)e^{i\omega k}

对于非周期信号,一般把整段信号视为一个周期,这种情况下的傅立叶变换称为离散时间傅立叶变换(Discrete Time Fourier Transform,DTFT)。类似地,DTFT 的公式为:

F(ω)=Tcn=k=f[k]eiωkF(\omega) =Tc_n=\sum_{k=-\infin}^{\infin}f[k]e^{-i\omega k}

逆变换为:

f[k]=12πππF(ω)eiωkdωf[k] =\frac{1}{2\pi}\int^{\pi}_{-\pi}F(\omega)e^{i\omega k} d\omega

📒结语

文章就到这里。限于水平,如果有纰漏之处,请不吝指出。

这篇博客稍微复习了一下傅立叶级数和傅立叶变换的数学推导,也补充了我以前没有注意到的相关原理,感觉日常使用的工具的原理清晰了许多。

0 条评论未登录用户
Ctrl or + Enter 评论
🌸 Run