如何避免FFT（matlab）计算傅里叶级数与傅里叶变换时存在的混叠（aliasing）和泄露（leakage）问题
文章目录
- 前言
- 1、 周期信号的傅里叶级数
- 1.1 周期信号的多种傅里叶表示形式
- 1.2 周期信号傅里叶级数系数的计算
- 2 非周期信号的傅里叶变换
- 2.1 非周期信号的傅里叶变换
- 2.2 傅里叶级数与傅里叶变换的关系
- 2.3 非周期信号的Parseval 能量等式
- 3 离散信号的傅里叶变换
- 3.1 离散傅里叶变换
- 3.2 快速傅里叶变换
- 3.3 快速傅里叶变换与傅里叶变换和傅里叶级数的关系
- 3.4 快速傅里叶变换的Parseval 能量等式
- 4信号时频变换中的失真问题
- 4.1 傅里叶级数中的失真问题
- 4.1.2 混叠
- 4.1.3 泄露
- 4.2 傅里叶变换中的失真问题
- 5 总结
前言
一个动力学信号的描述方式有两种，一种是时间域的描述，另外一种是频率域的描述。在时间域，以时间为横坐标，可以研究信号随时间的变化过程，如速度、加速度。在频域，以频率为横坐标，信号被不同频率的谐波描述，可以研究不同频率下信号幅值的分布。根据信号类型的不同（周期信号，非周期信号，离散信号等），时-频变换关系有所区别，涉及到傅里叶级数（Fourier Series，FS）、傅里叶变换（Fourier Transform，FT）、离散时间的傅里叶变换（discrete-time Fourier Transform，DTFT）等，而这些变换方法本身又有一定的关联。商业计算软件中，最常用的时-频处理方法是快速傅里叶变换方法（Fast Fourier Transform,FFT）。FFT与上述FS、FT和DTFT之间的关系，以及时-频变换中的信号失真问题是本文重点计算的内容。
本文的理论来源于Siu-Kui Au的《Operational Modal Analysis》，感兴趣的可以阅读原版推导。受作者水平所限，部分内容再加工过程中可能出现错误，欢迎批评指正。
1、 周期信号的傅里叶级数
1.1 周期信号的多种傅里叶表示形式
根据傅里叶理论，对于任何周期为T的周期信号x(t)x(t)x(t) 都可以表示为傅里叶级数（Fourier Series，FS）的形式
x(t)=a0+∑k=1∞akcosωˉkt+∑k=1∞bksinωˉktωˉk=2πkT(1)x\left(t\right)=a_0+\sum_{k=1}^{\infty}{a_k\cos{{\bar{\omega}}_kt}}+\sum_{k=1}^{\infty}{b_k\sin{{\bar{\omega}}_kt}}\ \ \ \ \ \ {\bar{\omega}}_k=\frac{2\pi k}{T} \tag{1} x(t)=a0+k=1∑∞akcosωˉkt+k=1∑∞bksinωˉkt ωˉk=T2πk(1)
傅里叶级数表示将周期信号分解为周期为 T/kT/kT/k 的无穷个简谐信号的和，当然也要考虑不随时间变化的常数项。对于特定的周期信号而言有限个简谐响应的和即可完成信号的模拟。
进一步可以将式（1）表示为幅值和相位的形式
x(t)=a0+∑k=1∞ak2+bk2cos(ωˉkt−ϕk)ϕk=arctanbkak(2)x\left(t\right)=a_0+\sum_{k=1}^{\infty}{\sqrt{a_k^2+b_k^2}\cos{\left({\bar{\omega}}_kt-\phi_k\right)}}\ \ \ \phi_k=\arctan{\frac{b_k}{a_k}} \tag{2} x(t)=a0+k=1∑∞ak2+bk2cos(ωˉkt−ϕk) ϕk=arctanakbk(2)
利用欧拉公式，可以将式（1）表示为复指数的形式
x(t)=c0+∑k=1∞ckeiωˉkt+∑k=1∞c−ke−iωˉkt(3)x\left(t\right)=c_0+\sum_{k=1}^{\infty}{c_ke^{i{\bar{\omega}}_kt}}+\sum_{k=1}^{\infty}{c_{-k}e^{-i{\bar{\omega}}_kt}} \tag{3} x(t)=c0+k=1∑∞ckeiωˉkt+k=1∑∞c−ke−iωˉkt(3)
进一步简写为
x(t)=∑k=−∞∞ckeiωˉktωˉk=2πkT(4)x\left(t\right)=\sum_{k=-\infty}^{\infty}{c_ke^{i{\bar{\omega}}_kt}}\ \ \ \ {\bar{\omega}}_k=\frac{2\pi k}{T} \tag{4} x(t)=k=−∞∑∞ckeiωˉkt ωˉk=T2πk(4)
复指数形式的傅里叶系数 {ck}k=−∞∞\left\{c_k\right\}_{k=-\infty}^\infty{ck}k=−∞∞也是复数，与式（1）的系数关系为
c0=a0ck=12(ak−ibk)c−k=12(ak+ibk)(5)c_0=a_0 \\\ c_k=\frac{1}{2}(a_k-ib_k) \\\ c_{-k}=\frac{1}{2}(a_k+ib_k) \tag{5} c0=a0 ck=21(ak−ibk) c−k=21(ak+ibk)(5)
1.2 周期信号傅里叶级数系数的计算
对式（1）在周期 t∈−T/2,T/2t\in{-T/2,T/2}t∈−T/2,T/2 内积分，可以计算常数 a0a_0a0 。分别对式（1）左右乘以 cosωˉkt\cos{{\bar{\omega}}_kt}cosωˉkt或者 sinωˉkt\sin{{\bar{\omega}}_kt}sinωˉkt ，并结合三角函数的正交性，可以计算系数 aka_kak 和 bkb_kbk ，计算结果为
a0=1T∫−T/2T/2x(t)dtak=2T∫−T/2T/2x(t)cosωˉktdtbk=2T∫−T/2T/2x(t)sinωˉktdt(6)a_0=\frac{1}{T}\int_{-T/2}^{T/2}{x(t)dt}\\\ a_k=\frac{2}{T}\int_{-T/2}^{T/2}{x(t)\cos{{\bar{\omega}}_kt}dt}\\\ b_k=\frac{2}{T}\int_{-T/2}^{T/2}{x(t)\sin{{\bar{\omega}}_kt}dt} \tag{6}a0=T1∫−T/2T/2x(t)dt ak=T2∫−T/2T/2x(t)cosωˉktdt bk=T2∫−T/2T/2x(t)sinωˉktdt(6)
通过式（5-6）可以计算得到
ck=1T∫−T/2T/2x(t)e−iωˉktdt(7)c_k=\frac{1}{T}\int_{-T/2}^{T/2}{x(t)e^{-i{\bar{\omega}}_kt}dt}\tag{7} ck=T1∫−T/2T/2x(t)e−iωˉktdt(7)
另外系数 c_k 具有共轭对称的性质即
c−k=cˉk(8)c_{-k}={\bar{c}}_k \tag{8} c−k=cˉk(8)
其中 xˉ{\bar{x}}xˉ 表示取共轭。
1.3 周期信号的Parseval能量等式
根据Parseval能量等式，信号的时域和频域应该具有相同能量。对于周期信号而言，在一个周期内其信号的能力是有限的，即∫−T/2T/2x2(t)dt<∞\int_{-T/2}^{T/2}{x^2(t)dt}<\infty∫−T/2T/2x2(t)dt<∞ 。其时-频域的能量等式为
1T∫−T/2T/2x2(t)dt=a02+12∑k=1∞(ak2+bk2)=∑k=−∞∞∣ck∣2(9)\frac{1}{T}\int_{-T/2}^{T/2}{x^2(t)dt}=a_0^2+\frac{1}{2}\sum_{k=1}^{\infty}{(a_k^2+b_k^2)}=\sum_{k=-\infty}^{\infty}\left|c_k\right|^2\tag{9} T1∫−T/2T/2x2(t)dt=a02+21k=1∑∞(ak2+bk2)=k=−∞∑∞∣ck∣2(9)
需要注意的是，上式表征的是时域平均能量与傅里叶系数的平方和的对等关系，可以理解为能量在单位时间与单位频点的平均。
2 非周期信号的傅里叶变换
2.1 非周期信号的傅里叶变换
对于非周期信号而言，并不能描述为傅里叶级数的形式，但是如果非周期信号的能量有限，依然可以写成简谐响应在连续频率区间的求和形式，被称为傅里叶变换（Fourier Transform，FT）。记x(t)x(t)x(t)为非周期信号，在 −∞<t<∞-\infty<t<\infty−∞<t<∞ 区间的能量有限，即 ∫−∞∞x2(t)dt<∞\int_{-\infty}^{\infty}{x^2(t)dt}<\infty∫−∞∞x2(t)dt<∞ ，则 x(t)x(t)x(t) 的傅里叶变换为
x(t)=12π∫−∞∞X(ω)eiωtdω(10)x\left(t\right)=\frac{1}{2\pi}\int_{-\infty}^{\infty}{X\left(\omega\right)e^{i\omega t}d\omega} \tag{10}x(t)=2π1∫−∞∞X(ω)eiωtdω(10)
X(ω)=∫−∞∞x(t)e−iωtdt(11)X\left(\omega\right)=\int_{-\infty}^{\infty}{x\left(t\right)e^{-i\omega t}dt} \tag{11}X(ω)=∫−∞∞x(t)e−iωtdt(11)
上式是傅里叶变换和傅里叶逆变换，与傅里叶级数的区别有两点：（1）傅里叶变换的频率是连续变量，傅里叶级数是离散变量；（2）傅里叶变换中有一个乘子 12π\frac{1}{2\pi}2π1 。
2.2 傅里叶级数与傅里叶变换的关系
从傅里叶级数出发，可以推导得到傅里叶变换。假设非周期的信号为 x(t)x(t)x(t) ，可以利用周期为 TTT 的周期信号 xp(t)x_p(t)xp(t) 近似表示，即 {x(t)=xp(t)t∈[−T/2,T/2]}\left\{\begin{matrix}x\left(t\right)=x_p(t)&t\in\left[-T/2,T/2\right]\\\end{matrix}\right\}{x(t)=xp(t)t∈[−T/2,T/2]} 。
当 T→∞T\rightarrow\inftyT→∞ ，则 xp(t)x_p(t)xp(t) 收敛于 x(t)x\left(t\right)x(t) ，对于给定的 TTT ，周期信号 xp(t)x_p(t)xp(t) 的复傅里叶系数为
ck=1T∫−T/2T/2xp(t)e−iωˉktdt=1T∫−T/2T/2x(t)e−iωˉktdt(12)c_k=\frac{1}{T}\int_{-T/2}^{T/2}{x_p(t)e^{-i{\bar{\omega}}_kt}dt}=\frac{1}{T}\int_{-T/2}^{T/2}{x\left(t\right)e^{-i{\bar{\omega}}_kt}dt} \tag{12}ck=T1∫−T/2T/2xp(t)e−iωˉktdt=T1∫−T/2T/2x(t)e−iωˉktdt(12)
根据Cauchy-Schwartz不等式有
∣ck∣2≤1T2∫−T2T2x2(t)dt×∫−T2T2∣e−iωˉkt∣2dt≤1T2∫−∞∞x2(t)dt×T=1T∫−∞∞x2(t)dt(13)\left|c_k\right|^2\le\frac{1}{T^2}\int_{-\frac{T}{2}}^{\frac{T}{2}}{x^2\left(t\right)dt}\times\int_{-\frac{T}{2}}^{\frac{T}{2}}{\left|e^{-i{\bar{\omega}}_kt}\right|^2dt}\le\frac{1}{T^2}\int_{-\infty}^{\infty}{x^2\left(t\right)dt}\times T=\frac{1}{T}\int_{-\infty}^{\infty}{x^2\left(t\right)dt} \tag{13}∣ck∣2≤T21∫−2T2Tx2(t)dt×∫−2T2T∣∣e−iωˉkt∣∣2dt≤T21∫−∞∞x2(t)dt×T=T1∫−∞∞x2(t)dt(13)
又因为 ∫−∞∞x2(t)dt<∞\int_{-\infty}^{\infty}{x^2\left(t\right)dt}<\infty∫−∞∞x2(t)dt<∞ ，则当 T→∞T\rightarrow\inftyT→∞ 时，∣ck∣→0\left|c_k\right|\rightarrow0∣ck∣→0 ，这也意味着对于非周期信号而言，不能使用傅里叶级数表示，因为乘子 1T\frac{1}{T}T1 的存在，使得傅里叶系数趋于无线小。
去掉傅里叶级数中的乘子1T\frac{1}{T}T1，并利用连续的频率替代离散频率，得到新的函数
Xp(ω)=∫−T2T2x(t)e−iωtdt(14)X_p\left(\omega\right)=\int_{-\frac{T}{2}}^{\frac{T}{2}}{x\left(t\right)e^{-i\omega t}dt} \tag{14} Xp(ω)=∫−2T2Tx(t)e−iωtdt(14)
则， ck=Xp(ωˉk)/Tc_k=X_p\left({\bar{\omega}}_k\right)/Tck=Xp(ωˉk)/T ，代入傅里叶级数公式有
xp(t)=∑k=−∞∞ckeiωˉkt=1T∑k=−∞∞Xp(ωˉk)eiωˉkt(15)x_p\left(t\right)=\sum_{k=-\infty}^{\infty}{c_ke^{i{\bar{\omega}}_kt}}=\frac{1}{T}\sum_{k=-\infty}^{\infty}{X_p\left({\bar{\omega}}_k\right)e^{i{\bar{\omega}}_kt}} \tag{15} xp(t)=k=−∞∑∞ckeiωˉkt=T1k=−∞∑∞Xp(ωˉk)eiωˉkt(15)
对于离散频点而言， Δω=ωˉk+1−ωˉk=2πT\Delta\omega={\bar{\omega}}_{k+1}-{\bar{\omega}}_k=\frac{2\pi}{T}Δω=ωˉk+1−ωˉk=T2π ，因此上式的1T\frac{1}{T}T1 可以改写为 1T=Δω2π\frac{1}{T}=\frac{\Delta\omega}{2\pi}T1=2πΔω。
又因为 T→∞T\rightarrow\inftyT→∞ ， xp(t)→x(t)x_p\left(t\right)\rightarrow x\left(t\right)xp(t)→x(t) ，则有
x(t)=limT→∞xp(t)=limT→∞∑k=−∞∞Xp(ωˉk)eiωˉktΔω2π(16)x\left(t\right)=\lim_{T\rightarrow\infty}{x_p(t)}=\lim_{T\rightarrow\infty}{\sum_{k=-\infty}^{\infty}{X_p\left({\bar{\omega}}_k\right)e^{i{\bar{\omega}}_kt}}\frac{\Delta\omega}{2\pi}} \tag{16}x(t)=T→∞limxp(t)=T→∞limk=−∞∑∞Xp(ωˉk)eiωˉkt2πΔω(16)
上式的右端可以改写为积分的形式为
x(t)=12π∫−∞∞X(ω)eiωtdω(17)x\left(t\right)=\frac{1}{2\pi}\int_{-\infty}^{\infty}{X\left(\omega\right)e^{i\omega t}d\omega} \tag{17}x(t)=2π1∫−∞∞X(ω)eiωtdω(17)
上述的推导过程中，我们交换了取极限和积分的顺序，在信号是有限能量的条件下，这是合理的。
2.3 非周期信号的Parseval 能量等式
对于非周期信号 x(t)x\left(t\right)x(t) 和其傅里叶变换形式 X(ω)X\left(\omega\right)X(ω) ，存在以下等式
∫−∞∞x2(t)dt=12π∫−∞∞∣X(ω)∣2dω(18)\int_{-\infty}^{\infty}{x^2(t)dt}=\frac{1}{2\pi}\int_{-\infty}^{\infty}{\left|X\left(\omega\right)\right|^2d\omega} \tag{18} ∫−∞∞x2(t)dt=2π1∫−∞∞∣X(ω)∣2dω(18)
3 离散信号的傅里叶变换
3.1 离散傅里叶变换
我们一般的采样信号都是离散时域信号，他在傅里叶级数和傅里叶变换中积分可以通过Riemann和的形式计算。记采样时间步长为 Δt\Delta tΔt，采样总数为 NNN ，信号 x(t)x(t)x(t) 的离散信号为 {xj=x(jΔt)}j=0N−1\left\{x_j=x\left(j\Delta t\right)\right\}_{j=0}^{N-1}{xj=x(jΔt)}j=0N−1 。则x(t)x(t)x(t)的傅里叶变化可以近似表示为
X(ω)=∫−∞∞x(t)eiωtdt≈X^(ω)=∑j=0N−1xje−iωjΔtΔt(19)X\left(\omega\right)=\int_{-\infty}^{\infty}{x\left(t\right)e^{i\omega t}dt}\approx\hat{X}\left(\omega\right)=\sum_{j=0}^{N-1}x_je^{-i\omega j\Delta t}\Delta t \tag{19} X(ω)=∫−∞∞x(t)eiωtdt≈X^(ω)=j=0∑N−1xje−iωjΔtΔt(19)
这里X^(ω)\hat{X}\left(\omega\right)X^(ω) 是 X(ω)X\left(\omega\right)X(ω) 的近似表达式，上式变换也被称为离散时间的傅里叶变换（discrete-time Fourier Transform，DTFT）， X^(ω)\hat{X}\left(\omega\right)X^(ω) 并不是对任何频点都能有效的估计，一般能够准确有效估计的频点为
ωk=2πkNΔtk=0,…,N−1(20)\omega_k=\frac{2\pi k}{N\Delta t}\ \ k=0,\ldots,N-1 \tag{20}ωk=NΔt2πk k=0,…,N−1(20)
因此 X^(ω)\hat{X}\left(\omega\right)X^(ω) 也被写为 X^(ωk)\hat{X}\left(\omega_k\right)X^(ωk) ，注意在傅里叶级数中，我们用到的频点符号是 ωˉk=2πkT{\bar{\omega}}_k=\frac{2\pi k}{T}ωˉk=T2πk 。
3.2 快速傅里叶变换
离散信号在 ωk\omega_kωk 频点处的DTFT可以通过快速傅里叶变换（Fast Fourier Transform，FFT）计算。记离散信号{xj}j=0N−1的FFT序列为{yk}j=0N−1\left\{x_j\right\}_{j=0}^{N-1}的FFT序列为\left\{y_k\right\}_{j=0}^{N-1}{xj}j=0N−1的FFT序列为{yk}j=0N−1， ，它们的对应关系可以定义为
yk=∑j=0N−1xje−i2πjk/N(21)y_k=\sum_{j=0}^{N-1}{x_je^{-i2\pi jk/N}} \tag{21}yk=j=0∑N−1xje−i2πjk/N(21)
{yk}j=0N−1\left\{y_k\right\}_{j=0}^{N-1}{yk}j=0N−1 的逆FFT序列定义为 {zj}j=0N−1\left\{z_j\right\}_{j=0}^{N-1}{zj}j=0N−1 ，定义为
zj=1N∑k=0N−1ykei2πjk/N(22)z_j=\frac{1}{N}\sum_{k=0}^{N-1}{y_ke^{i2\pi jk/N}} \tag{22}zj=N1k=0∑N−1ykei2πjk/N(22)
这里需要注意的是 zjz_jzj 是复数值，而 xjx_jxj 是实数值。另外 yky_kyk 也被称为离散傅里叶变换（discrete Fourier Transform，DFT），由于其计算方式通常使用FFT，因此本文不区分这两者的微小区别。
{yk}j=0N−1\left\{y_k\right\}_{j=0}^{N-1}{yk}j=0N−1 也有镜像共轭的特点，即
yN−k=ykˉ(23)y_{N-k}=\bar{y_k} \tag{23} yN−k=ykˉ(23)
MATLAB中的函数fft得到计算结果为yky_kyk，且满足共轭特征。
3.3 快速傅里叶变换与傅里叶变换和傅里叶级数的关系
时域信号的傅里叶变换为 X(ω)X\left(\omega\right)X(ω) ，离散信号可以近似表示为 X^(ω)\hat{X}\left(\omega\right)X^(ω) ，在特定的有效频点 ωk\omega_kωk 处近似表示为 X^(ωk)\hat{X}\left(\omega_k\right)X^(ωk) ，简记为 X^k{\hat{X}}_kX^k。将 ωk=2πkNΔt\omega_k=\frac{2\pi k}{N\Delta t}ωk=NΔt2πk 和 yky_kyk 代入，有
X(ωk)≈X^k=ykΔt(24)X\left(\omega_k\right)\approx{\hat{X}}_k=y_k\Delta t \tag{24} X(ωk)≈X^k=ykΔt(24)
对于傅里叶级数，假设 {xj}j=0N−1\left\{x_j\right\}_{j=0}^{N-1}{xj}j=0N−1 是周期为 TTT 的周期信号的离散采样信号， NΔt=TN\Delta t=TNΔt=T ，则 ωk=2πkNΔt\omega_k=\frac{2\pi k}{N\Delta t}ωk=NΔt2πk
与 ωˉk=2πkT{\bar{\omega}}_k=\frac{2\pi k}{T}ωˉk=T2πk 对应。同样利用Riemann求和计算傅里叶级数的积分，可以得到
ck=1T∫−T/2T/2x(t)e−iωˉktdt≈c^k=1NΔt∑j=0N−1xje−i2πjkN×Δt=ykN(25)c_k=\frac{1}{T}\int_{-T/2}^{T/2}{x\left(t\right)e^{-i{\bar{\omega}}_kt}dt}\approx{\hat{c}}_k=\frac{1}{N\Delta t}\sum_{j=0}^{N-1}{x_je^{-\frac{i2\pi jk}{N}}}\times\Delta t=\frac{y_k}{N} \tag{25} ck=T1∫−T/2T/2x(t)e−iωˉktdt≈c^k=NΔt1j=0∑N−1xje−Ni2πjk×Δt=Nyk(25)
为了方便于在实际应用中（当周期未知时或者希望使用多个周期估计数据的傅里叶级数来平均噪声），这里并不是严格要求 NΔt=TN\Delta t=TNΔt=T 。
3.4 快速傅里叶变换的Parseval 能量等式
对于 {xj}j=0N−1和{yk}j=0N−1\left\{x_j\right\}_{j=0}^{N-1} 和 \left\{y_k\right\}_{j=0}^{N-1}{xj}j=0N−1和{yk}j=0N−1 ，有如下等式
∑j=0N−1xj2=1N∑j=0N−1∣yk∣2(26)\sum_{j=0}^{N-1}x_j^2=\frac{1}{N}\sum_{j=0}^{N-1}\left|y_k\right|^2 \tag{26}j=0∑N−1xj2=N1j=0∑N−1∣yk∣2(26)
4信号时频变换中的失真问题
傅里叶级数和傅里叶变换的离散时间近似导致了信号的特征性误差。只要来自于两个方面：混叠（aliasing）和泄露（leakage）。混叠指的是信号中高出Nyquist frequency（ 1/2Δt1/2\Delta t1/2Δt ）的高频部分被错误的混在了低频信号分析中。一个简单的介绍，人的眼睛接受视觉信号的频率是有限的（就像采样频率一样），当汽车速度较慢时，我们能清晰的看到车轮辐条的转动速度，但是当汽车速度较快时，我们会发现我们看到汽车辐条的转动速度与真实速度并不符合，有时候甚至是静止的，这就是混叠现象。混叠现象在采样完成之后是无法避免的，只能事先过滤掉，信号中超过Nyquist frequency的谐波成分；泄露的发生的主要原因是信号在测量时域内没有整数倍周期特性，通常可以通过增加持续时间来抑制泄露的发生。本节分别介绍傅里叶级数与傅里叶转化中的失真问题。
4.1 傅里叶级数中的失真问题
假设周期信号为 x(t)x(t)x(t) ，其复傅里叶系数为 {ck}k=−∞∞\left\{c_k\right\}_{k=-\infty}^\infty{ck}k=−∞∞ 。周期信号对应的离散时间序列为 {xj=x(jΔt)}j=0N−1\left\{x_j=x\left(j\Delta t\right)\right\}_{j=0}^{N-1}{xj=x(jΔt)}j=0N−1 ，离散时间信号对应的FFT系数为 {yk}k=0N−1\left\{y_k\right\}_{k=0}^{N-1}{yk}k=0N−1 ， ckc_kck 的离散时间近似为 c^k=yk/N{\hat{c}}_k=y_k/Nc^k=yk/N 。
4.1.2 混叠
当信号 x(t)x(t)x(t) 中含有超过Nyquist frequency的谐波成分时，会发生混叠现象。假设数据持续时间等于信号周期，则 c^k{\hat{c}}_kc^k 为
c^k=ck+∑m=1∞(cmN+k+cmN−kˉ)(27){\hat{c}}_k=c_k+\sum_{m=1}^{\infty}\left(c_{mN+k}+\bar{c_{mN-k}}\right) \tag{27}c^k=ck+m=1∑∞(cmN+k+cmN−kˉ)(27)
其中 c^k{\hat{c}}_kc^k 不仅包含频点 fkf_kfk 的贡献，还包含了频点 fs±fk，2fs±fkf_s\pm f_k，{2f}_s\pm f_kfs±fk，2fs±fk，…， 的贡献。 fs=1/Δtf_s=1/\Delta tfs=1/Δt 为采样频率， fk=k/NΔtf_k=k/N\Delta tfk=k/NΔt 为FFT变换的主要频点。混叠发生在采样时间是信号周期的整数倍条件下，当采样时间不是信号周期的整数倍时，还会发生信号泄露现象。
算例1：
考虑一个周期信号 x(t)=2cos2πftx\left(t\right)=2\cos{2\pi ft}x(t)=2cos2πft ，傅里叶级数的系数 ak,bk{a_k,b_k}ak,bk 只有 a1=2a_1=2a1=2 ，其余都为0。复傅里叶系数 {ck}\left\{c_k\right\}{ck} 中，只有 c1=1,c−1=1c_1=1 , c_{-1}=1c1=1,c−1=1 。 假设离散数据 {xj}j=0N−1\left\{x_j\right\}_{j=0}^{N-1}{xj}j=0N−1 的采样时间间隔 Δt=0.1s\Delta t=0.1sΔt=0.1s ，采样数目N=10N=10N=10 ，采样时长 t=NΔt=1st=N\Delta t=1st=NΔt=1s 。通过FFT方法计算复傅里叶级数的系数为 c^k=yk/N{\hat{c}}_k=y_k/Nc^k=yk/N 。
| 图1 离散周期信号的混叠 |
|---|
图1中，超过Nyquist frequency即5Hz的数据是其小于5Hz数据的共轭，这里为了更好的说明问题也一并画出，实际应用中可以省去。当信号的周期不超过Nyquist frequency即5Hz时，可以得到准确的傅里叶变换结果。但是当频率超过5 Hz时，会发生混叠现象，高频成分会影响低频频点的准确性。 f=6Hzf=6\ \rm{Hz}f=6 Hz 时，则在4Hz处会混叠在fs−fk=10−4=6Hzfs-f_k=10-4=6\ \rm Hzfs−fk=10−4=6 Hz 。
clc,clear,close all
delta_t = 0.1;
fs = 1/delta_t;
N = 10;
t = N*delta_t;
t_array = (0:N-1)*delta_t;
f_array = (0:N-1)/t;f = 1;
x_array = 2*cos(2*pi*f*t_array);
c_k = fft(x_array)/N;subplot(1,2,1)
title('fs = 1 Hz')
t_actual = 0:0.001:1;
x_actual = 2*cos(2*pi*f*t_actual);
plot(t_actual,x_actual,'--')
hold on
plot(t_array,x_array,'bo','MarkerSize',6,...'MarkerEdgeColor','b',...'MarkerFaceColor','b')
xlabel('t(s)','FontSize',15)
title(['f = ',num2str(f),'Hz'],'FontSize',15)
subplot(1,2,2)
stem(f_array,c_k,'filled','MarkerSize',6,...'MarkerEdgeColor','b',...'MarkerFaceColor','b')
xlim([0 10])
ylim([0 1.2])
xlabel('f_k','FontSize',15)
4.1.3 泄露
当信号采样的时间 NΔtN\Delta tNΔt 不是信号周期 TTT 的整倍数时，傅里叶级数的FFT近似表示会发生泄露现象。这种条件下，傅里叶级数的频点 ωˉr=2πrT{\bar{\omega}}_r=\frac{2\pi r}{T}ωˉr=T2πr 与FFT的频点 ωk=2πkNΔt\omega_k=\frac{2\pi k}{N\Delta t}ωk=NΔt2πk 不匹配，傅里叶级数的谐响应会泄露到其他的FFT频点处。FFT近似的复傅里叶级数系数可以表述为复傅里叶级数系数与核函数的卷积形式
KaTeX parse error: \tag works only in display equations
其中
K1(ω)=1N∑j=0N−1e−iωjΔt=sin(NωΔt/2)Nsin(ωΔt/2)e−i(N−1)ωΔt/2(29)K_1\left(\omega\right)=\frac{1}{N}\sum_{j=0}^{N-1}e^{-i\omega j\Delta t}=\frac{\sin{(N\omega\Delta t/2)}}{N\sin{(\omega\Delta t/2)}}e^{-i(N-1)\omega\Delta t/2} \tag{29}K1(ω)=N1j=0∑N−1e−iωjΔt=Nsin(ωΔt/2)sin(NωΔt/2)e−i(N−1)ωΔt/2(29)
令 u=ωΔt2πu=\frac{\omega\Delta t}{2\pi}u=2πωΔt ，则 K1(ω)K_1\left(\omega\right)K1(ω) 的无量纲表述形式为
DN(u)=sin(Nπu)Nsin(πu)(31)D_N\left(u\right)=\frac{\sin{(N\pi u)}}{N\sin{(\pi u)}} \tag{31}DN(u)=Nsin(πu)sin(Nπu)(31)
DN(u)D_N\left(u\right)DN(u) 被称为狄利克雷核函数。
上式说明了 c^k{\hat{c}}_kc^k 包含了所有傅里叶级数频点的贡献。来自频点 ωˉr{\bar{\omega}}_rωˉr 的贡献不直接来自 crc_rcr ，需要一个衰减因子 K1(ωk−ωˉr)K_1(\omega_k-{\bar{\omega}}_r)K1(ωk−ωˉr) 的参与，而衰减的强度取决于 ωˉr{\bar{\omega}}_rωˉr 到 ωk\omega_kωk 的距离。
| 图2 狄利克雷核函数 |
|---|
图2 是 ∣DN(u)∣\left|D_N\left(u\right)\right|∣DN(u)∣ 示意图。周期为1，N=20N=20N=20，主瓣为 [−12,12][-\frac{1}{2},\frac{1}{2}][−21,21] ，当 u∈[−12,12]u\in[-\frac{1}{2},\frac{1}{2}]u∈[−21,21] 时， ω∈[−πΔt,πΔt]\omega\in[-\frac{\pi}{\Delta t},\frac{\pi}{\Delta t}]ω∈[−Δtπ,Δtπ] ，即 $frac{\pi}{\Delta t}$ 为Nyquist频点。从图中可以看出，当 u=±1N,±2Nu=\pm\frac{1}{N},\pm\frac{2}{N}u=±N1,±N2 等FFT主要频点处， DN(u)D_N\left(u\right)DN(u) 为0。
由于卷积效应， c^k{\hat{c}}_kc^k 包含的谐波是周期信号 x(t)x(t)x(t) 在以下频点的数据：
（1）FFT主要频点 fk=kNΔtHzf_k=\frac{k}{N\Delta t} \rm Hzfk=NΔtkHz 。
（2） fkf_kfk 附近的 O(fsN)HzO(\frac{f_s}{N}) \rm HzO(Nfs)Hz （泄露）
（3） fs±fk,2fs±fkf_s\pm f_k,2f_s\pm f_kfs±fk,2fs±fk ,…附近的 O(fsN)HzO(\frac{f_s}{N}) \rm HzO(Nfs)Hz （混叠+泄露）
算例2
按照算例1的信号，只是采样点数目调整为 N=12N=12N=12 ，相应的采样时长为 NΔt=1.2sN\Delta t=1.2\ \rm sNΔt=1.2 s 。计算结果如图3所示，对于频率为1、4、6和11Hz的信号而言，采样时长与信号周期比分别为1.2、2.4、7.2和13.2，都不是整数倍周期，因此会导致信号泄露的发生，另外当频率为6Hz和11Hz时，还有信号混叠现象。
| 图3 离散周期信号的泄露现象 |
|---|
4.2 傅里叶变换中的失真问题
受Nyquist的限制，混叠和泄露也存在于傅里叶变换的FFT近似表达中，机理与傅里叶级数的失真类似。
将时域离散信号 xj=12π∫−∞∞X(ω′)eiω′jΔtdω′x_j=\frac{1}{2\pi}\int_{-\infty}^{\infty}{X\left(\omega^\prime\right)e^{i\omega^\prime j\Delta t}d\omega^\prime}xj=2π1∫−∞∞X(ω′)eiω′jΔtdω′ 代入 X^(ω)\hat{X}(\omega)X^(ω) 的近似公式中，得到
X^(ω)=∑j=0N−1[12π∫−∞∞X(ω′)eiω′jΔtdω′]e−iω′jΔtΔt=12π∫−∞∞X(ω′)NΔt×1N∑j=0N−1e−i(ω−ω′)jΔtdω′=12π∫−∞∞X(ω′)K2(ω−ω′)dω′(32)\hat{X}\left(\omega\right)=\sum_{j=0}^{N-1}{\left[\frac{1}{2\pi}\int_{-\infty}^{\infty}{X\left(\omega^\prime\right)e^{i\omega^\prime j\Delta t}d\omega^\prime}\right]e^{-i\omega^\prime j\Delta t}\Delta t} \\\ =\frac{1}{2\pi}\int_{-\infty}^{\infty}{X\left(\omega^\prime\right)N\Delta t\times\frac{1}{N}\sum_{j=0}^{N-1}e^{-i{(\omega-\omega}^\prime)j\Delta t}d\omega^\prime}\\\ =\frac{1}{2\pi}\int_{-\infty}^{\infty}{X\left(\omega^\prime\right)K_2({\omega-\omega}^\prime)d\omega^\prime} \tag{32} X^(ω)=j=0∑N−1[2π1∫−∞∞X(ω′)eiω′jΔtdω′]e−iω′jΔtΔt =2π1∫−∞∞X(ω′)NΔt×N1j=0∑N−1e−i(ω−ω′)jΔtdω′ =2π1∫−∞∞X(ω′)K2(ω−ω′)dω′(32)
其中
K2(ω)=NΔtK1(ω)=NΔtDN(u)e−iπ(N−1)uu=ωΔt2π(33)K_2\left(\omega\right)=N\Delta tK_1\left(\omega\right)=N\Delta tD_N\left(u\right)e^{-i\pi\left(N-1\right)u}\ \ \ \ u=\frac{\omega\Delta t}{2\pi} \tag{33}K2(ω)=NΔtK1(ω)=NΔtDN(u)e−iπ(N−1)u u=2πωΔt(33)
从上式可以看出， X^(ω)\hat{X}(\omega)X^(ω) 也是 X(ω′)X\left(\omega^\prime\right)X(ω′) 与核函数的卷积的形式。根据该式也可以解释傅里叶变化中的混叠和泄露现象。对于混叠而言，如图2所示 ∣DN(u)∣\left|D_N\left(u\right)\right|∣DN(u)∣ 在 u=0,±1,±2,…u=0,\pm1,\pm2,\ldotsu=0,±1,±2,… 附件存在局部最大值1，对应的函数 ∣K2(ω−ω′)∣\left|K_2({\omega-\omega}^\prime)\right|∣K2(ω−ω′)∣ 在 ω′=ω±2πmfsm=0,1,2,…\omega^\prime=\omega\pm2\pi mf_s\ \ \ \ m=0,1,2,\ldotsω′=ω±2πmfs m=0,1,2,… 也存在局部最大值。那么 X^(ω)\hat{X}(\omega)X^(ω) 的频点贡献来源不仅是 ω\omegaω ，而是 ω±2πmfs\omega\pm2\pi mf_sω±2πmfs 。对于泄露， X^(ω)\hat{X}(\omega)X^(ω) 在 ω\omegaω 频点的贡献来源于原信号在 ω′\omega^\primeω′ 频点的傅里叶变化值 X(ω′)X\left(\omega^\prime\right)X(ω′) ，并且这个贡献被算子 K2(ω−ω′)K_2({\omega-\omega}^\prime)K2(ω−ω′) 衰减，衰减强度根据两频点的距离 (ω−ω′)({\omega-\omega}^\prime)(ω−ω′) 决定。傅里叶级数中泄露的原因是采样时间不是信号周期的整倍数导致的，但是在傅里叶变换中有所不同。傅里叶变换的信号是非周期信号，频点是连续的，这就意味着无论采样时间无法满足信号周期的整倍数，因此傅里叶变换中的泄露问题也就无视采样时长了，然而，随着采样时间的增加，泄露的影响会越来越小，因为当 ∣ω−ω′∣≫2π/NΔt\left|{\omega-\omega}^\prime\right|\gg2\pi/N\Delta t∣ω−ω′∣≫2π/NΔt 时， ∣K2(ω−ω′)∣\left|K_2({\omega-\omega}^\prime)\right|∣K2(ω−ω′)∣ 的值会很小。
5 总结
本文前述中涉及到四种‘傅里叶’：傅里叶级数、傅里叶变换、离散时间的傅里叶变换以及快速傅里叶变换。针对不同的信号形式，需要采用不同的傅里叶形式，其中最基本的是针对周期信号的傅里叶级数形式，其基本思想是将周期信号写成无穷个离散频率的简谐响应之和，傅里叶级数有三种种表述方法，分别为实傅里叶级数系数、傅里叶级数相位和幅值和复傅里叶级数系数。如果信号是非周期信号，在满足能量有限的前提，其依然可以写成各频点谐响应求和的形式，但是这里的频点是连续的，意味着每个频点对非周期信号的贡献都是有限的，这种变换形式被称为傅里叶数变换，非周期信号可以视为周期无限长的周期信号，因此傅里叶变换的推导来源于傅里叶级数。然而实际科学应用中，时间信号以时域离散信号为主，因此利用黎曼求和公式处理傅里叶变换中涉及到的积分问题，得到离散时域信号的离散近似傅里叶变换，也被称为离散时间傅里叶变换，而离散时间变换中最高效的计算方法为快速傅里叶变换方法，该方法也是各商业软件中最常用的方法。利用快速傅里叶变换方法可以近似计算傅里叶级数和傅里叶变换的频域信号，但是可计算的有效频点并不是任意的，和采样时间间隔和采样时长有关系。
离散信号的傅里叶级数和傅里叶变换都是近似结果，因此与实际结果会存在一定的误差，产生误差的原因来混叠和泄露。然而这两个因素的影响又和采样频率和采样时长息息相关。当信号的频率成分有有超过Nyquist频率（采样频率的一半）的信号源时，超过的部分不能被有效计算，且会混叠在小于Nyquist频率的某些频点中，造成整个傅里叶变换的失真，这就是混叠现象，这种现象可以通过加大采样频率和滤掉高频信号源消除。对于泄露问题，周期信号的傅里叶级数和非周期信号的傅里叶变换有一点区别。周期信号的泄露是因为采样时长与周期信号的周期不是整倍数关系，导致了傅里叶级数频点与信号本身的频点不同步，可以通过调整采样周期来缓解周期信号的信号泄露问题。然而非周期信号的傅里叶数变换中，频点是连续的，信号泄露是不可避免的，但是可以通过增加采样时长来有效减小信号泄露对信号失真的影响。
如何避免FFT（matlab）计算傅里叶级数与傅里叶变换时存在的混叠（aliasing）和泄露（leakage）问题相关推荐
- dft计算傅里叶级数系数_DFT（离散傅里叶变换）与FFT（快速傅里叶变换）初识
一. 简介 离散傅里叶变换(Discrete Fourier Transform, DFT)是数字信号处理最重要的基石之一,也是对信号进行分析和处理时最常用的工具之一.在200多年前法国数学家.物理学 ...
- 基于MATLAB去理解掌握傅里叶级数和傅里叶变换
1.周期信号的傅里叶级数 f(t)=f(t+T) F0=1/T为基波频率 满足狄利赫里条件则周期信号可以展开为三角函数的线性组合 (1) 在一个周期内,函数f(t)为连续或只含有有限个第一类间断点: ...
- 用MATLAB计算序列的离散傅里叶变换
用MATLAB计算序列的离散傅里叶变换 MATLAB提供了用快速算法计算离散傅里叶变换的函数fft,其调用格式为: Xk = fft(xn, N) 其中,调用参数xn为时域序列向量,N为离散傅里叶变换 ...
- 卷积、傅里叶级数、傅里叶变换、快速傅里叶变换、pytorch中的fft，rfft
卷积: 连续形式: 离散形式: '卷' : 翻转 和 滑动 '积' : 积分 翻转:g(t) - > g(-t) 滑动:g(-t) - > g(n-t) 平移n个单位 举个例 ...
- python计算无穷级数求和常用公式_傅里叶变换（二） 从傅里叶级数到傅里叶变换...
在上一部分当中,得到了利用三角函数表示周期函数的方法,但是对于非周期函数就...凉了.所以有什么办法吗?没办法(划掉).这时候我们就需要拿出来我们的黑科技--傅里叶变换. 一.傅里叶级数的推广 当然这 ...
- 《数字图像处理》-（3）-1从傅里叶级数到傅里叶变换详细推导以及傅里叶图像的性质
1 傅里叶级数到傅里叶变换公式推导 1.1傅里叶级数 傅里叶级数:周期信号都可以分解为有限或无限个正弦波或余弦波的叠加,且这些波的频率都是原始信号频率的整数倍.用傅里叶级数或变换表示的函数特征完全可以 ...
- Matlab：二维傅里叶变换
Matlab:二维傅里叶变换 二维傅里叶变换 二维衍射模式 fft2 函数将二维数据变换为频率空间.例如,您可以变换二维光学掩膜以揭示其衍射模式. 二维傅里叶变换 以下公式定义 m×n 矩阵 X 的离 ...
- matlab对一组数据傅里叶变换,matlab快速傅里叶变换（三个matlab程序介绍）
描述 一种积分变换,它来源于函数的傅里叶积分表示.积分 (1) 称为ƒ 的傅里叶积分.周期函数在一定条件下可以展成傅里叶级数,而在(-∞,∞)上定义的非周期函数ƒ,显然不能用三角级数来表示.但是J.- ...
- 时频分析方法总结：傅里叶级数及傅里叶变换、STFT 、小波变换、Wigner-Ville 分布
前言: 一.傅里叶变换的机理 一个能量无限的正弦信号和源信号乘积并求和得到某个频率下的系数,随着频率的增加,正弦信号改变,再次求得系数,依次构成了频谱图 傅里叶级数及傅里叶变换 https://blo ...
最新文章
- 北大数学天才许晨阳，回国效力6年后，为什么又去了美国任教？
- 生成对抗网络学习笔记4----GAN（Generative Adversarial Nets）的实现
- 【ACL 2020】腾讯AI Lab解读三大前沿方向及入选的20篇论文
- android 定义date对象,如何从Date对象设置Android Chronometer基准时间？
- DPDK — 网卡初始化流程（Intel 82599 ixgbe 网卡驱动示例）
- 提升购物体验，跨境电商如何做企业管理？
- python高效编程15个利器_你不知道的18个Python高效编程技巧
- 大数据WEB阶段Spring框架（二）简化配置的操作
- python-random种子
- SSRS 动态设置分组依据及行组个数
- bind()的实现(持续更新中)
- linux tc 对本机网卡限速
- TextView跑马灯和editText抢占焦点，键盘弹不出来问题解决
- 自制1寸照片及打印排版
- 九大Java性能优化工具
- vscode 显示/设置隐藏文件夹
- Pytorch每日一练——预测泰坦尼克号船上的生存乘客
- 测试-嵌入式图床外链
- 嵌入式Linux驱动编程复习资料
- 万里长征第一步——Hello World