二维傅里叶变换(2D Fourier Transforms)可视化

By Long Luo

前言

前几天看了一个介绍 核磁共振成像(MRI) 原理的视频:
三维演示磁共振成像(MRI)原理
,加上自己也曾做过 22 次核磁共振,一直很好奇其原理,所以花了点时间弄懂了核磁共振成像原理。

在这篇文章我们不详细介绍核磁共振成像原理,只关注其成像部分,也就是根据获取到的空间频率来还原图像,也就是二维傅里叶逆变换(2D Inverse Fourier Transforms\textit{2D Inverse Fourier Transforms})

关于傅里叶变换及傅里叶逆变换,可以参考之前写过的一篇文章:傅里叶变换 。傅里叶变换本质上就是换基,推荐这篇文章:我理解的傅里叶变换

一维傅里叶变换(1D Fourier Transform)

对于非周期函数,我们把它的周期看作无穷大。基频 ω0=2πT\omega_0=\frac{2\pi}{T} 则趋近于无穷小,又因为基频相当于周期函数的傅里叶级数中两个相邻频率的差值 (n+1)ω0nω0(n+1)\omega_0-n\omega_0,我们可以把它记作 Δω\Delta\omega 或者微分 dωd\omeganω0n\omega_0 则相当于连续变量 ω\omega。这样就得到了针对非周期函数的频谱函数如下:

cn=Δω2π+f(t)ejwtdtc_{n}=\frac{\Delta\omega}{2\pi}\int_{-\infty}^{+\infty} f(t)e^{-jwt}dt

我们会发现这里的 cnc_n 是趋于 00 的。

将它代入 f(t)=n=n=+cnejnω0tf(t)=\sum_{n=-\infty}^{n=+\infty}c_ne^{jn\omega_0t}

得到:

f(t)=+(Δω2π+f(t)ejwtdt)ejωt=+12π(+f(t)ejωtdt)ejωtdωf(t)=\sum_{-\infty}^{+\infty}(\frac{\Delta\omega}{2\pi}\int_{-\infty}^{+\infty} f(t)e^{-jwt}dt)e^{j\omega t}=\int_{-\infty}^{+\infty}\frac{1}{2\pi}(\int_{-\infty}^{+\infty}f(t)e^{-j\omega t}dt)e^{j\omega t}d\omega

则非周期函数的傅里叶变换定义为

F(ω)=+f(t)ejωtdtF(\omega)=\int_{-\infty}^{+\infty}f(t)e^{-j\omega t}dt

我们可以发现 cn=dω2πF(ω)c_n=\frac{d\omega}{2\pi}F(\omega) ,即我们选取了信号幅值在频域中的分布密度 F(ω)F(\omega) 来表示傅里叶变换,而不是相应频率的信号幅值大小 cnc_n 。因为如果选择后者,那傅里叶变换的函数值就都是无穷小了,这显然对我们没有任何帮助。

我们一般也用频率 ff 来进行傅里叶变换:

F(f)=+f(t)ej2πftdtF(f)=\int_{-\infty}^{+\infty}f(t)e^{-j2\pi f t}dt

在数学上 f(t)f(t) 大多是在时域 tt 上是连续的,但是由于计算机采集信号在时域中是离散的,所以实际应用中的 f(t)f(t) 都是其经采样处理后得到的 fs(t)f_s(t)

同时,计算机也只能计算出有限个频率上对应的幅值密度,即 ff 也是离散的。

DFT(Discrete Fourier Transform)\textit{DFT(Discrete Fourier Transform)} 也就是 ttff 都是离散的傅里叶变换。

一维离散傅里叶变换

采样(Sampling)

首先引入冲激函数(也叫 Dirac\textit{Dirac} 函数)的概念。

t0t \neq 0δ(t)=0\delta(t)=0+δ(t)dt=1\int_{-\infty}^{+\infty}\delta(t)dt=1

根据它的定义,我们可知 +δ(tt0)f(t)dt=f(t0)\int_{-\infty}^{+\infty}\delta(t-t_0)f(t)dt=f(t_0) 。这是 Dirac\textit{Dirac} 函数的重要性质,容易看出,它可以筛选出 f(t)f(t)t0t_0 处的函数值,即起到采样的作用。

但是 Dirac\textit{Dirac} 函数一次只能选取一个函数值,所以我们将很多个采样点不同的 Dirac\textit{Dirac} 函数叠加起来,就可以实现时域上的采样了。这样叠加的函数被称为梳状函数,表达式为:

δs(t)=n=δ(tnTs)\delta_s(t)=\sum_{n=-\infty}^{\infty}\delta(t-nT_s) ,其中 TsT_s 为采样周期。

将时域上的连续信号与它相乘,即可得到 fs=n=f(t)δ(tnTs)f_s=\sum_{n=-\infty}^{\infty}f(t)\delta(t-nT_s)

时域离散化计算

对于采样得到的 xs(t)=n=x(t)δ(tnTs)x_s(t)=\sum_{n=-\infty}^{\infty}x(t)\delta(t-nT_s) 进行傅里叶变换。

根据傅里叶变换的定义

X(ω)=+x(t)ejωtdtX(\omega)=\int_{-\infty}^{+\infty}x(t)e^{-j\omega t}dt

Xs(ω)=(n=x(t)δ(tnTs))ejωtdtX_s(\omega)=\int_{-\infty}^{\infty}(\sum_{n=-\infty}^{\infty}x(t)\delta(t-nT_s))e^{-j\omega t}dt

交换积分与求和顺序,得到

Xs(ω)=n=x(t)δ(tnTs)ejwtdtX_s(\omega)=\sum_{n=-\infty}^{\infty}\int_{-\infty}^{\infty}x(t)\delta(t-nT_s)e^{-jwt}dt

根据 Dirac\textit{Dirac} 函数的筛选特性, Xs(ω)=n=x(nTs)ejwnTsX_s(\omega)=\sum_{n=-\infty}^{\infty}x(nT_s)e^{-jwnT_s}

这样就完成了我们的时域离散化计算。

频域离散化计算

时域离散化得到的结果 Xs(ω)X_s(\omega) 在频域上仍是连续的,而计算机只能求取有限个 ω\omega 对应的频谱密度。此外, Xs(ω)X_s(\omega) 中的时域采样次数 NN 为无穷大,实际应用中显然不会进行无穷多次时域采样。

我们首先解决 NN 为无穷大的问题。对于连续信号 x(t)x(t) 进行 NN 次(NN 为有限值)采样,采样周期为 TsT_s 。然后对采样得到的信号进行时域上的周期延拓(延拓至正负无穷),这样我们就得到了一个周期为 T0=NTsT_0=NT_s 的函数。对于周期函数而言,它的频谱密度函数是离散化的,这样我们就成功把频域也进行了离散化。具体计算方法如下:

在一个周期内,离散信号的表达式为 xs(t)=n=0N1x(t)δ(tnTs)x_s(t)=\sum_{n=0}^{N-1}x(t)\delta(t-nT_s)

离散信号的傅里叶级数:

X(kω0)=1T00T(n=0N1x(t)δ(tnTs))ejkw0tdtX(k\omega_0)=\frac{1}{T_0}\int_{0}^{T}(\sum_{n=0}^{N-1}x(t)\delta(t-nT_s))e^{-jkw_0t}dt ,其中 ω0=2πT0\omega_0=\frac{2\pi}{T_0}

交换积分和求和次序,

X(kω0)=1T0n=0N10Tx(t)δ(tnTs)ejkω0tdtX(k\omega_0)=\frac{1}{T_0}\sum_{n=0}^{N-1}\int_{0}^{T}x(t)\delta(t-nT_s)e^{-jk\omega_0t}dt

根据狄拉克函数的筛选特性, X(kω0)=1T0n=0N1x(nTs)ejkω0nTsX(k\omega_0)=\frac{1}{T_0}\sum_{n=0}^{N-1}x(nTs)e^{-jk\omega_0nT_s}

X(kω0)=1NTsn=0N1x(nTs)ej2πNTsknTs=1NTsn=0N1x(nTs)ej2πNknX(k\omega_0)=\frac{1}{NT_s}\sum_{n=0}^{N-1}x(nTs)e^{-j\frac{2\pi}{NT_s}knT_s}=\frac{1}{NT_s}\sum_{n=0}^{N-1}x(nTs)e^{-j\frac{2\pi}{N}kn}

为更加简明,令 X[k]=X(kω0)TsX[k]=X(k\omega_0)T_sx[n]=x(nTs)x[n]=x(nT_s),则

X[k]=1Nn=0N1x[n]ej2πNknX[k]=\frac{1}{N}\sum_{n=0}^{N-1}x[n]e^{-j\frac{2\pi}{N}kn}

其中:NN 为采样总点数,nn 为当前采样点,x[n]x[n] 在时域 nn 的信号值,kk 为当前频率值(从 00 Hz 到 N1N-1 Hz),X[k]X[k] DFT\textit{DFT} 的结果(是一个复数,包含幅度和相位信息)。

一维傅里叶变换可视化

square wave

FFT and IFFT

440Hz Music

FFT 440Hz

Mixed Music

FFT Mixed Music

二维连续傅里叶变换(2D Fourier Transforms)

二维连续傅里叶变换

F(u,v)=f(x,y)ej2π(ux+vy)dxdyF(u,v) = \int_{-\infty }^{\infty}\int_{-\infty }^{\infty} f(x,y)e^{-j2\pi(ux+vy)}dxdy

二维连续傅里叶逆变换

f(x,y)=14π2F(u,v)ej2π(ux+vy)dudvf(x,y) = \frac{1}{4\pi^2} \int_{-\infty }^{\infty}\int_{-\infty }^{\infty} F(u,v)e^{j2\pi(ux+vy)}dudv

二维离散傅里叶变换

二维离散傅里叶变换

F(u,v)=x=0M1y=0N1f(x,y)ej2π(uxM+vyN)F(u,v)=\sum_{x=0}^{M-1} \sum_{y=0}^{N-1}f(x,y)e^{-j2\pi(\frac{ux}{M}+\frac{vy}{N})}

二维离散傅里叶逆变换

f(x,y)=1MNx=0M1y=0N1F(u,v)ej2π(uxM+vyN)f(x,y)=\frac{1}{MN}\sum_{x=0}^{M-1} \sum_{y=0}^{N-1}F(u,v)e^{j2\pi(\frac{ux}{M}+\frac{vy}{N})}

傅里叶变换特征参数

F(u,v)=R(u,v)+jI(u,v)F(u,v)=R(u,v)+jI(u,v)

频谱/幅度谱

F(u,v)=R2(u,v)+I2(u,v)|F(u,v)|=\sqrt{R^2(u,v)+I^2(u,v)}

P(u,v)=R2(u,v)+I2(u,v)P(u,v)=R^2(u,v)+I^2(u,v)

相位谱

ψ(u,v)=arctanI(u,v)R(u,v)\psi (u,v)=arctan\frac{I(u,v)}{R(u,v)}

参考资料