跳至主要內容

离散傅里叶变换

大约 10 分钟

离散傅里叶变换

前提结论

傅里叶变换的对称性

对于傅里叶变换对 f(t)f(t) 与 F(f)F(f) (注意此处采用的频域自变量单位为频率, 对于角频率还需要而外修正)

F[f(t)]=F(f) \mathscr{F}[f(t)]=F(f)

假设有信号具有与频谱相同的函数 F(t)F(t)
则该信号的傅里叶变换满足

F[F(t)]=f(−f) \mathscr{F}[F(t)]=f(-f)

傅里叶级数与变换的联系

无论是傅里叶级数还是变换, 信号的频域特性均不变
因此对于周期信号 f(t)f(t) 的傅里叶级数与变换间满足关系 (注意 f0f_0 为[基频], 单位为频率)

Cn=F(nf0)df C_n=F(nf_0)\mathrm{d}f

根据冲激函数的性质, 可得周期信号的傅里叶变换为

F(f)=∑n=−∞∞Cnδ(f−nf0) F(f)=\sum_{n=-\infty}^{\infty}C_n\delta(f-nf_0)

傅里叶变换的离散性

由傅里叶级数与变换的联系得, 周期信号 f(t)f(t) 的傅里叶变换 F(f)F(f) 必定是离散的, 且间隔为 f0f_0, 即信号的频率
根据对称性还可以推出, 对于离散信号 fs(t)f_s(t) 的傅里叶变换 Fs(f)F_s(f) 必定是周期性的, 且周期为 1Ts\frac{1}{T_s}, 即离散信号间隔 TsT_s 的倒数

脉冲信号的傅里叶变换

对于周期为 TsT_s 的脉冲信号

δs(t)=∑n=−∞∞δ(t−nTs) \delta_s(t)=\sum_{n=-\infty}^{\infty}\delta(t-nT_s)

其有复数傅里叶级数

Cn=1Ts∫−Ts2Ts2δs(t)e−j2πnt/Tsdt=1Ts C_n=\frac{1}{T_s}\int_{-\frac{T_s}{2}}^{\frac{T_s}{2}}\delta_s(t)e^{-j2\pi nt/T_s}\mathrm{d}t=\frac{1}{T_s}

因此 Cn=1TsC_n=\frac{1}{T_s} 为一个与 nn 无关的常数, 根据傅里叶级数与变换的联系, 可得其傅里叶变依然为一个脉冲函数

F(f)=1Ts∑n=−∞∞δ(f−nfs)=1Tsδfs(t) F(f)=\frac{1}{T_s}\sum_{n=-\infty}^{\infty}\delta(f-nf_s)=\frac{1}{T_s}\delta_{fs}(t)

周期延拓

对于信号 f(t)f(t), 根据卷积的特性, 将其与间隔为 T0T_0 的脉冲函数 δ0(t)\delta_0(t) 求卷积得到信号 f0(t)f_0(t) 必定为以 T0T_0 周期
特别的, 当 f(t)f(t) 为宽度为 T0T_0 的有限信号, f0(t)f_0(t) 为 f(t)f(t) 在空间上的重复

(f∗δ0)(t)=∫−∞∞δs(τ)f(t−τ)dτ=∫−∞∞∑n=−∞∞f(t−τ)δ(τ−nT0)dτ=∑n=−∞∞f(t−nT0) \begin{split}(f*\delta_0)(t)&=\int_{-\infty}^{\infty}\delta_s(\tau)f(t-\tau)\mathrm{d}\tau\\ &=\int_{-\infty}^{\infty}\sum_{n=-\infty}^{\infty}f(t-\tau)\delta(\tau-nT_0)\mathrm{d}\tau\\ &=\sum_{n=-\infty}^{\infty}f(t-nT_0) \end{split}

约定

  • 以下推导中采用 sinc⁡(t)=sin⁡tt\operatorname{sinc}(t)=\frac{\sin t}{t}

信号采样

定义

  • 时域上周期无限长的信号 f(t)f(t)
  • 采样间隔为 TsT_s 的信号采样脉冲 δs(t)=∑n=−∞∞δ(t−nTs)\delta_s(t)=\sum_{n=-\infty}^{\infty}\delta(t-nT_s)
  • 截取信号的窗函数 u0(t)=h(t+Ts2)−h(t−Ts2−T0)u_0(t)=h(t+\frac{T_s}{2})-h(t-\frac{T_s}{2}-T_0)

通过采样脉冲对信号进行采样, 得到间隔为 TsT_s 的离散信号

fs(t)=f(t)δs(t) f_s(t)=f(t)\delta_s(t)

仅采样 t=−Ts/2∼T0−Ts/2t=-T_s/2\sim T_0-T_s/2 间共 NN 个信号, 注意采样间隔 TsT_s 与截取长度 T0T_0 间满足 Ts=NT0T_s=NT_0
相当于使用一个矩形窗函数 u0(t)u_0(t) 截取信号

f^(t)=fs(t)u0(t)=∑n=0N−1f(t)δ(t−nTs) \hat{f}(t)=f_s(t)u_0(t)=\sum_{n=0}^{N-1}f(t)\delta(t-nT_s)

此时 f^(t)\hat{f}(t) 便为一个可通过计算机处理的理想有限离散信号

频域采样

对于有限离散信号 f^(t)\hat{f}(t), 其傅里叶变换 F^(f)\hat{F}(f) 必定为连续的周期函数, 依然不利于计算机处理
因此也需要对信号的傅里叶变换做采样

根据周期延拓的性质, 傅里叶变换最适合以信号宽度的倒数 f0f_0 为间隔采样有

F~(f)=F^(f)δ0(f) \tilde{F}(f)=\hat{F}(f)\delta_0(f)

由傅里叶变换的卷积特性以及脉冲函数的傅里叶变换可得, 此频域采样也等价于将有限信号 f^(t)\hat{f}(t) 在空间上延拓, 每个周期 −Ts/2∼T0+Ts/2-T_s/2\sim T_0+T_s/2 内即 f^(t)\hat{f}(t)

f~(t)=f^(t)∗F−1[δ0(f)]=T0∑r=−∞∞f^(t−rT0)=T0∑r=−∞∞∑n=0N−1f(t−rT0)δ(t−rT0−nTs)=T0∑r=−∞∞∑n=0N−1f(nTs)δ(t−rT0−nTs) \begin{split}\tilde{f}(t)&=\hat{f}(t)*\mathscr{F}^{-1}[\delta_0(f)]\\ &=T_0\sum_{r=-\infty}^{\infty}\hat{f}(t-rT_0)\\ &=T_0\sum_{r=-\infty}^{\infty}\sum_{n=0}^{N-1}f(t-rT_0)\delta(t-rT_0-nT_s)\\ &=T_0\sum_{r=-\infty}^{\infty}\sum_{n=0}^{N-1}f(nT_s)\delta(t-rT_0-nT_s) \end{split}

此时 f~(t)\tilde{f}(t), 有周期 T0T_0, 可求出其傅里叶级数 (积分时直接取 r=0,f^(t)r=0, \hat{f}(t) 相同的周期)

Ck=1T0∫−Ts/2T0+Ts/2T0∑n=0N−1f(nTs)δ(t−nTs)e−j2πkt/T0dt=∑n=0N−1f(nTs)e−j2πknTs/T0=∑n=0N−1f(nTs)e−j2πkn/N \begin{split}C_k&=\frac{1}{T_0}\int_{-T_s/2}^{T_0+T_s/2}T_0\sum_{n=0}^{N-1}f(nT_s)\delta(t-nT_s)e^{-j2\pi kt/T_0}\mathrm{d}t\\ &=\sum_{n=0}^{N-1}f(nT_s)e^{-j2\pi knT_s/T_0}\\ &=\sum_{n=0}^{N-1}f(nT_s)e^{-j2\pi kn/N} \end{split}

根据傅里叶级数与变换的联系, 可由此得到其傅里叶变换

F~(f)=∑k=−∞∞∑n=0N−1f(nTs)e−j2πkn/Nδ(f−kf0)(a) \tilde{F}(f)=\sum_{k=-\infty}^{\infty}\sum_{n=0}^{N-1}f(nT_s)e^{-j2\pi kn/N}\delta(f-kf_0)\tag{a}

注意, 根据 e−j2π(k+N)n/N=e−j2π(n/N+1)=1⋅e−j2πn/Ne^{-j2\pi (k+N)n/N}=e^{-j2\pi (n/N+1)}=1\cdot e^{-j2\pi n/N} (或傅里叶变换的离散性) 可得, F~(f)\tilde{F}(f) 依然为一个周期函数, 周期为 fs=Nf0f_s=Nf_0

因此使用离散序列 F~[k]\tilde{F}[k] 表示为

F~[k]=F~(kf0)df=∑n=0N−1f(nTs)e−j2πkn/N(b) \tilde{F}[k]=\tilde{F}(kf_0)\mathrm{d}f=\sum_{n=0}^{N-1}f(nT_s)e^{-j2\pi kn/N}\tag{b}

逆变换

通过以上推导, 来自 f(t)f(t) 的周期离散信号 f~(t)\tilde{f}(t) 的傅里叶变换 F[f~(t)]=F~(f)\mathscr{F}[\tilde{f}(t)]=\tilde{F}(f) 已经求得, 现求其逆变换 F−1[F~(f)]=f~(t)\mathscr{F}^{-1}[\tilde{F}(f)]=\tilde{f}(t)

不难发现, F~(f)\tilde{F}(f) 与 f~(t)\tilde{f}(t) 均为离散周期函数
因此可以尝试将 F~(f)\tilde{F}(f) 的周期分离, 使其与 f~(t)\tilde{f}(t) 具有相同的形式
之后根据傅里叶变换的对称性有 F[F~(t)]=f~(−f)\mathscr{F}[\tilde{F}(t)]=\tilde{f}(-f)

对于序列 F~[k]=F~(kf0)\tilde{F}[k]=\tilde{F}(kf_0) 具有周期 NN, 因此可令 r=k+rN,k=0∼N−1r=k+rN, k=0\sim N-1
现对 F~(t)\tilde{F}(t) 进行变形

F~(f)=∑k=−∞∞∑n=0N−1f(nTs)e−j2πkn/Nδ(f−kf0)=∑r=−∞∞∑k=0N−1∑n=0N−1f(nTs)e−j2π(k+rN)n/Nδ(f−kf0−rNf0)=fs∑r=−∞∞∑k=0N−1[1fs∑n=0N−1f(nTs)e−j2πkn/N]δ(f−kf0−rfs) \begin{split}\tilde{F}(f)&=\sum_{k=-\infty}^{\infty}\sum_{n=0}^{N-1}f(nT_s)e^{-j2\pi kn/N}\delta(f-kf_0)\\ &=\sum_{r=-\infty}^{\infty}\sum_{k=0}^{N-1}\sum_{n=0}^{N-1}f(nT_s)e^{-j2\pi (k+rN)n/N}\delta(f-kf_0-rNf_0)\\ &=f_s\sum_{r=-\infty}^{\infty}\sum_{k=0}^{N-1}\Bigg[\frac{1}{f_s}\sum_{n=0}^{N-1}f(nT_s)e^{-j2\pi kn/N}\Bigg]\delta(f-kf_0-rf_s) \end{split}

不难发现, 除 1fs∑n=0N−1f(nTs)e−j2πkn/N\frac{1}{f_s}\sum_{n=0}^{N-1}f(nT_s)e^{-j2\pi kn/N} 部分外, 此时的 F~(t)\tilde{F}(t) 具有与 f~(t)\tilde{f}(t) 完全相同的形式

根据离散序列 F~[k]\tilde{F}[k] 的定义可得, 此部分即

Ts∑n=0N−1f(nTs)e−j2πkn/N=TsF~[k] T_s\sum_{n=0}^{N-1}f(nT_s)e^{-j2\pi kn/N}=T_s\tilde{F}[k]

使用 TsF~[n]T_s\tilde{F}[n] 代替 f(nTs)f(nT_s), TsT_s 代替 f0f_0, 带入式 a, 然后令 k=−k,t=−tk=-k, t=-t 有

f~(−t)=∑k=−∞∞∑n=0N−1TsF~[n]e−j2πkn/Nδ(t−kTs)f~(t)=Ts∑k=−∞∞∑n=0N−1F~[n]e−j2π(−k)n/Nδ(−t+kTs)f~(t)=Ts∑k=−∞∞∑n=0N−1F~[n]ej2πkn/Nδ(t−kTs) \begin{split}\tilde{f}(-t)&=\sum_{k=-\infty}^{\infty}\sum_{n=0}^{N-1}T_s\tilde{F}[n]e^{-j2\pi kn/N}\delta(t-kT_s)\\ \tilde{f}(t)&=T_s\sum_{k=-\infty}^{\infty}\sum_{n=0}^{N-1}\tilde{F}[n]e^{-j2\pi (-k)n/N}\delta(-t+kT_s)\\ \tilde{f}(t)&=T_s\sum_{k=-\infty}^{\infty}\sum_{n=0}^{N-1}\tilde{F}[n]e^{j2\pi kn/N}\delta(t-kT_s) \end{split}

同样的, 使用离散序列 f~[k]\tilde{f}[k] 表示为

f~[k]=f~(kTs)dt=Ts∑n=0N−1F~[n]ej2πkn/N \tilde{f}[k]=\tilde{f}(kT_s)\mathrm{d}t=T_s\sum_{n=0}^{N-1}\tilde{F}[n]e^{j2\pi kn/N}

离散傅里叶级数 DFS

将式 b带入 f~[k]\tilde{f}[k] 中的 F~[n]\tilde{F}[n] 有 (或根据频域采样步骤中 f(t)f(t) 与 f~(t)\tilde{f}(t) 的关系可得到相同结论)

f~[k]=Ts∑n=0N−1∑m=0N−1f(mTs)e−j2πnm/Nej2πkn/N=Ts∑n=0N−1∑m=0N−1f(mTs)ej2πn(k−m)/N=Ts∑n=0N−1f(kTs)=NTsf(kTs) \begin{split}\tilde{f}[k]&=T_s\sum_{n=0}^{N-1}\sum_{m=0}^{N-1}f(mT_s)e^{-j2\pi nm/N}e^{j2\pi kn/N}\\ &=T_s\sum_{n=0}^{N-1}\sum_{m=0}^{N-1}f(mT_s)e^{j2\pi n(k-m)/N}\\ &=T_s\sum_{n=0}^{N-1}f(kT_s)\\ &=NT_sf(kT_s) \end{split}

综上所述, 有离散傅里叶级数 (DFSDFS) 变换对

{F~[k]=1NTs∑n=0N−1f~[n]e−j2πkn/N,k∈Zf~[k]=Ts∑n=0N−1F~[n]ej2πkn/N,k∈Z \begin{cases} \tilde{F}[k]=\frac{1}{NT_s}\sum_{n=0}^{N-1}\tilde{f}[n]e^{-j2\pi kn/N},&k\in Z\\ \\ \tilde{f}[k]=T_s\sum_{n=0}^{N-1}\tilde{F}[n]e^{j2\pi kn/N},&k\in Z \end{cases}

序列 F~[k]\tilde{F}[k] 与 f~[k]\tilde{f}[k] 均为无限长序列, 且具有周期 NN, 分别有分布间隔 f0=fs/Nf_0=f_s/N 与 TsT_s

离散傅里叶变换 DFT

在实际使用中, 通常会直接使用采集到的离散数据序列 f[k]f[k], 其定义即为

f[k]=f(kTs),k∈[0,N) f[k]=f(kT_s), k\in[0,N)

由于计算机中的序列不可能无限长, 因此仅截取序列 F~[k]\tilde{F}[k] 一个周期 NN 内的值
定义 f[k]f[k] 的离散傅里叶变换为序列 F[k]F[k] 满足

F[k]=F~[k],k∈[0,N) F[k]=\tilde{F}[k], k\in[0,N)

规定 Ts=1T_s=1, 使用 f[k]f[k] 代替离散傅里叶级数中的 f~[k]=NTsf[k]\tilde{f}[k]=NT_sf[k]
因此有离散傅里叶变换 (DFTDFT) 变换对

{F[k]=∑n=0N−1f[n]e−j2πkn/N,k∈[0,N)f[k]=1N∑n=0N−1F[n]ej2πkn/N,k∈[0,N) \begin{cases} F[k]=\sum_{n=0}^{N-1}f[n]e^{-j2\pi kn/N},&k\in [0,N)\\ \\ f[k]=\frac{1}{N}\sum_{n=0}^{N-1}F[n]e^{j2\pi kn/N},&k\in [0,N) \end{cases}

信号复原

复原信号时, 可令 kTs=tkT_s=t, 得到复原信号

f′(t)=1N∑n=0N−1F[n]ej2πtn/T0 f'(t)=\frac{1}{N}\sum_{n=0}^{N-1}F[n]e^{j2\pi tn/T_0}

由于 F[k]F[k] 也经过 f=−f0/2∼fs−f0/2f=-f_0/2\sim f_s-f_0/2 的截取
因此复原得到的实际信号为 (对信号取绝对值)

f′(t)=∣f~(t)∗[fssinc⁡(πfst)e−j2πfct/2]∣=f~(t)∗fssinc⁡(πfst) f'(t)=|\tilde{f}(t)*[f_s\operatorname{sinc}(\pi f_s t)e^{-j2\pi f_ct/2}]|=\tilde{f}(t)*f_s\operatorname{sinc}(\pi f_s t)

其中 fc=(fs−f0)/2f_c=(f_s-f_0)/2 可根据傅里叶变换的移频性自行推导

根据卷积特性可得 f′(t)=fs∑k=−∞∞f[k]s(t−kTs)f'(t)=f_s\sum_{k=-\infty}^{\infty}f[k]s(t-kT_s)
其中 s(t)=sinc⁡(πfst)s(t)=\operatorname{sinc}(\pi f_s t), 且 s(nTs)=0(n≠0)/=1(n=0)s(\frac{n}{T_s})=0(n\neq 0)/=1(n=0)

因此 k∈[0,N)k\in [0,N) 时必然有 f′(kT)=fsf(kT)f'(kT)=f_sf(kT) 但其余部分无法保证, 但只要采样点 NN 足够多, 就能复原出原信号
根据以上推导可得, F[k]F[k] 复原得到的是 fsf(kT)f_s f(kT), 因此 F[k]F[k] 反应的是 fsf(t)f_s f(t) 的频域特性
因此根据傅里叶变换的线性性, 对于 FFTFFT 的结果 F[t]F[t], 还需要取 F′[k]=F[k]/fsF'[k]=F[k]/f_s 才能得到接近实际傅里叶变换 F(t)F(t) 的结果

DFT 与 DFS 的关系

以下部分为个人猜想, 有待证明

通过离散傅里叶级数 DFSDFS 得到的结果为离散的幅值谱, 近似认为

F[k]≈fs⋅∫(k−12)f0(k+12)f0F(f)df F[k]\approx f_s\cdot \int_{(k-\frac{1}{2})f_0}^{(k+\frac{1}{2})f_0}F(f)\mathrm{d}f

对于 F[k]F[k] 与 F~[k]\tilde{F}[k] 可发现当输入序列相同时满足

F~[k]=f0⋅F[k] \tilde{F}[k]=f_0\cdot F[k]

通过离散傅里叶变换 DFTDFT 得到的结果为离散的幅值密度谱, 近似认为 (待证明)

F[k]≈T0⋅fs⋅∫(k−12)f0(k+12)f0F(f)df F[k]\approx T_0\cdot f_s\cdot \int_{(k-\frac{1}{2})f_0}^{(k+\frac{1}{2})f_0}F(f)\mathrm{d}f

因此 DFTDFT 结果的最大分辨率为 1T0\frac{1}{T_0}

以下推论经实验为正确的, 但成因不确定, 也可能由能量泄露导致

由上可知对于如 sin⁡,cos⁡\sin,\cos 等幅值密度谱中有冲击函数的信号
假设其在 (k−12)f0∼(k+12)f0{(k-\frac{1}{2})f_0}\sim{(k+\frac{1}{2})f_0} 的傅里叶变换结果为 F(f)=Aδ(kf0)F(f)=A\delta(kf_0)
此时其 DFTDFT 的结果为 F[k]=T0fsAF[k]=T_0 f_s A

DFT 特点总结

  1. 采样长度 T0T_0, 采样间隔 TsT_s, 采样点数 N=T0/TsN=T_0/T_s
  2. 首先使用窗函数截取连续信号 f(t)f(t) 在 t=−Ts/2∼T0−Ts/2t=-T_s/2\sim T_0-T_s/2 的部分, 此时将导致能量泄露
  3. 然后使用间隔为 TsT_s 的脉冲采样, 此时频域以 1Ts\frac{1}{T_s} 为周期循环, 当信号频域有 f>fs/2f>f_s/2 的部分, 将与其他周期混淆, 即香农采样定理
  4. 之后以间隔为 1T0\frac{1}{T_0} 的脉冲在频域采样, 因此 DFTDFT 结果的最大分辨率为 1T0\frac{1}{T_0}, 此时时域的采样信号将以 T0T_0 为周期循环
  5. 对于 DFTDFT 得到的结果, 还需要除以 fsf_s 才能与接近实际的傅里叶变换结果对应

实际使用

实际中多使用快速傅里叶变换 (FFTFFT), 其算法与离散傅里叶变换相同, 但是优化了运算过程
具体体现为

  1. 将截取区间修改为 k∈(−N2,N2]k\in (-\frac{N}{2},\frac{N}{2}]
  2. 当数据点不为 2n2^n 时, 将自动补 00

对于 numpy.fft 中的 f = fft.fft(x) 函数
参数 xx 为实数数组, 即采集到的离散数据点
返回值 f
f 为复数数组, 即离散傅里叶变换结果
索引 [0:N/2] 保存了正数部分的幅值 F[0]∼F[N/2]F[0]\sim F[N/2]
索引 [N/2+1:] 保存了负数部分的幅值 F[−N/2]∼F[−1]F[-N/2]\sim F[-1]

根据信号复原可得, 变换结果 F[k]F[k] 对应的是频率为 f=kT0=kNTsf=\frac{k}{T_0}=\frac{k}{NT_s} 下的幅值

误差分析

f′(t)f'(t) 的形式也表明, 经过离散傅里叶变换后, 仅能保留频率为 f=nT0f=\frac{n}{T_0} 的成分

注意

该笔记尚未完成