By Long Luo
前言
前几天看了一个介绍 核磁共振成像(MRI) 原理的视频:
三维演示磁共振成像(MRI)原理 ,加上自己也曾做过 2 次核磁共振,一直很好奇其原理,所以花了点时间弄懂了核磁共振成像原理。
在这篇文章我们不详细介绍核磁共振成像原理,只关注其成像部分,也就是根据获取到的空间频率来还原图像,也就是二维傅里叶逆变换(2D Inverse Fourier Transforms)。
关于傅里叶变换及傅里叶逆变换,可以参考之前写过的一篇文章:傅里叶变换 。傅里叶变换本质上就是换基,推荐这篇文章:我理解的傅里叶变换 。
对于非周期函数,我们把它的周期看作无穷大。基频 ω0=T2π 则趋近于无穷小,又因为基频相当于周期函数的傅里叶级数中两个相邻频率的差值 (n+1)ω0−nω0,我们可以把它记作 Δω 或者微分 dω。nω0 则相当于连续变量 ω。这样就得到了针对非周期函数的频谱函数如下:
cn=2πΔω∫−∞+∞f(t)e−jwtdt
我们会发现这里的 cn 是趋于 0 的。
将它代入 f(t)=∑n=−∞n=+∞cnejnω0t
得到:
f(t)=−∞∑+∞(2πΔω∫−∞+∞f(t)e−jwtdt)ejωt=∫−∞+∞2π1(∫−∞+∞f(t)e−jωtdt)ejωtdω
则非周期函数的傅里叶变换定义为
F(ω)=∫−∞+∞f(t)e−jωtdt
我们可以发现 cn=2πdωF(ω) ,即我们选取了信号幅值在频域中的分布密度 F(ω) 来表示傅里叶变换,而不是相应频率的信号幅值大小 cn 。因为如果选择后者,那傅里叶变换的函数值就都是无穷小了,这显然对我们没有任何帮助。
我们一般也用频率 f 来进行傅里叶变换:
F(f)=∫−∞+∞f(t)e−j2πftdt
在数学上 f(t) 大多是在时域 t 上是连续的,但是由于计算机采集信号在时域中是离散的,所以实际应用中的 f(t) 都是其经采样处理后得到的 fs(t) 。
同时,计算机也只能计算出有限个频率上对应的幅值密度,即 f 也是离散的。
DFT(Discrete Fourier Transform) 也就是 t 和 f 都是离散的傅里叶变换。
一维离散傅里叶变换
采样(Sampling)
首先引入冲激函数(也叫 Dirac 函数)的概念。
当 t=0 时 δ(t)=0; ∫−∞+∞δ(t)dt=1。
根据它的定义,我们可知 ∫−∞+∞δ(t−t0)f(t)dt=f(t0) 。这是 Dirac 函数的重要性质,容易看出,它可以筛选出 f(t) 在 t0 处的函数值,即起到采样的作用。
但是 Dirac 函数一次只能选取一个函数值,所以我们将很多个采样点不同的 Dirac 函数叠加起来,就可以实现时域上的采样了。这样叠加的函数被称为梳状函数,表达式为:
δs(t)=∑n=−∞∞δ(t−nTs) ,其中 Ts 为采样周期。
将时域上的连续信号与它相乘,即可得到 fs=∑n=−∞∞f(t)δ(t−nTs) 。
时域离散化计算
对于采样得到的 xs(t)=∑n=−∞∞x(t)δ(t−nTs) 进行傅里叶变换。
根据傅里叶变换的定义
X(ω)=∫−∞+∞x(t)e−jωtdt
则 Xs(ω)=∫−∞∞(∑n=−∞∞x(t)δ(t−nTs))e−jωtdt
交换积分与求和顺序,得到
Xs(ω)=n=−∞∑∞∫−∞∞x(t)δ(t−nTs)e−jwtdt
根据 Dirac 函数的筛选特性, Xs(ω)=∑n=−∞∞x(nTs)e−jwnTs
这样就完成了我们的时域离散化计算。
频域离散化计算
时域离散化得到的结果 Xs(ω) 在频域上仍是连续的,而计算机只能求取有限个 ω 对应的频谱密度。此外, Xs(ω) 中的时域采样次数 N 为无穷大,实际应用中显然不会进行无穷多次时域采样。
我们首先解决 N 为无穷大的问题。对于连续信号 x(t) 进行 N 次(N 为有限值)采样,采样周期为 Ts 。然后对采样得到的信号进行时域上的周期延拓(延拓至正负无穷),这样我们就得到了一个周期为 T0=NTs 的函数。对于周期函数而言,它的频谱密度函数是离散化的,这样我们就成功把频域也进行了离散化。具体计算方法如下:
在一个周期内,离散信号的表达式为 xs(t)=∑n=0N−1x(t)δ(t−nTs)
离散信号的傅里叶级数:
X(kω0)=T01∫0T(∑n=0N−1x(t)δ(t−nTs))e−jkw0tdt ,其中 ω0=T02π。
交换积分和求和次序,
X(kω0)=T01∑n=0N−1∫0Tx(t)δ(t−nTs)e−jkω0tdt
根据狄拉克函数的筛选特性, X(kω0)=T01∑n=0N−1x(nTs)e−jkω0nTs
即
X(kω0)=NTs1n=0∑N−1x(nTs)e−jNTs2πknTs=NTs1n=0∑N−1x(nTs)e−jN2πkn
为更加简明,令 X[k]=X(kω0)Ts , x[n]=x(nTs),则
X[k]=N1n=0∑N−1x[n]e−jN2πkn
其中:N 为采样总点数,n 为当前采样点,x[n] 在时域 n 的信号值,k 为当前频率值(从 0 Hz 到 N−1 Hz),X[k] DFT 的结果(是一个复数,包含幅度和相位信息)。
一维傅里叶变换可视化






二维连续傅里叶变换
F(u,v)=∫−∞∞∫−∞∞f(x,y)e−j2π(ux+vy)dxdy
二维连续傅里叶逆变换
f(x,y)=4π21∫−∞∞∫−∞∞F(u,v)ej2π(ux+vy)dudv
二维离散傅里叶变换
二维离散傅里叶变换
F(u,v)=x=0∑M−1y=0∑N−1f(x,y)e−j2π(Mux+Nvy)
二维离散傅里叶逆变换
f(x,y)=MN1x=0∑M−1y=0∑N−1F(u,v)ej2π(Mux+Nvy)
傅里叶变换特征参数
F(u,v)=R(u,v)+jI(u,v)
频谱/幅度谱
∣F(u,v)∣=R2(u,v)+I2(u,v)
P(u,v)=R2(u,v)+I2(u,v)
相位谱
ψ(u,v)=arctanR(u,v)I(u,v)
参考资料