DFT/IDFT 离散傅里叶转换/逆转换及代码实现 笔记

2020/4/16 傅里叶变换

离散傅里叶变换(DFT),是傅里叶变换在时域和频域上都呈现离散的形式,将时域信号的采样变换为在离散时间傅里叶变换(DTFT)频域的采样。

# 定义

# DFT (离散傅里叶转换)

XF(k)=n=0N1x(n)e(j2πN)kn,k=0,1,2,,N1X^F(k) = \sum_{n=0}^{N-1} x(n) e^{(\frac{-j2\pi}{N})kn}, \qquad k=0,1,2,\cdots,\ N-1 \\

x(n)x(n) 是一个均匀采样序列,即均匀时间间隔,n=0,1,,N1n=0,1,\cdots,N-1XF(k)X^F(k) 为第 KK 个 DFT 系数,k=0,1,,N1k=0,1,\cdots,N-1j=1j=\sqrt{-1}

# IDFT (离散傅里叶转换)

x(n)=1Nk=0N1XF(k)WNkn,k=0,1,2,,N1 为、N 点采样数据(2.1b)x(n) = \frac{1}{N}\sum_{k=0}^{N-1} X^F(k) W_N^{-kn}, \qquad k=0,1,2,\cdots,\ N-1\ 为、N\ 点采样数据 \tag{2.1b}

(WNkn)=WNkn=exp(j2πNkn),e±jθ=cos(θ)±jsin(θ)(W_N^{kn})^* = W_N^{-kn} = \exp(\frac{j2\pi}{N}kn), \qquad e^{\pm j\theta} = \cos(\theta) \pm j\sin(\theta)

其中, * 表示复共轭运算。DFT 的变换对可以使用下式表示:

x(n)XF(k)x(n)\Leftrightarrow X^F(k)

式 (2.1b) 中的归一化因子 1N\dfrac{1}{N} 可以平均分配在 DFT 和 IDFT 中(称为归一化 DFT)。也可以全部放在 DFT 运算中。

# 归一化 DFT

{XF(k)=1Nn=0N1x(n)WNknk=0,1,,N1x(n)=1Nk=0N1XF(k)WNknn=0,1,,N1\left\{ \begin{aligned} X^{\mathrm{F}}(k)=\frac{1}{\sqrt{N}} \sum_{n=0}^{N-1} x(n) W_{N}^{k n} \qquad k=0,1, \cdots, N-1 \\ x(n)=\frac{1}{\sqrt{N}} \sum_{k=0}^{N-1} X^{\mathrm{F}}(k) W_{N}^{-k n} \quad n=0,1, \cdots, N-1 \end{aligned} \right.

{XF(k)=1Nn=0N1x(n)exp(j2πknN)k=0,1,,N1x(n)=k=0N1XF(k)exp(j2πknN)n=0,1,,N1\left\{ \begin{aligned} X^{\mathrm{F}}(k)&=\frac{1}{N} \sum_{n=0}^{N-1} x(n) \exp \left(\frac{-\mathrm{j} 2 \pi k n}{N}\right) \quad &k=0,1, \cdots, N-1 \\ x(n)&=\sum_{k=0}^{N-1} X^{\mathrm{F}}(k) \exp \left(\frac{\mathrm{j} 2 \pi k n}{N}\right) \qquad &n=0,1, \cdots, N-1 \end{aligned} \right.

利用欧拉公式:e±jθ=cos(θ)±jsin(θ)e^{\pm j\theta} = \cos(\theta) \pm j\sin(\theta) ,有

{XF(k)=n=0N1x(n)[cos(2πknN)jsin(2πknN)]k=0,1,,N1x(n)=1Nk=0N1XF(k)[cos(2πknN)+jsin(2πknN)]n=0,1,,N1(2.2)\left\{ \begin{aligned} X^{\mathrm{F}}(k)&=\sum_{n=0}^{N-1} x(n) [\cos(\frac{2\pi kn}{N})-j\sin(\frac{2\pi kn}{N})] \quad &k=0,1, \cdots, N-1 \\ x(n)&=\frac{1}{N} \sum_{k=0}^{N-1} X^{\mathrm{F}}(k) [\cos(\frac{2\pi kn}{N})+j\sin(\frac{2\pi kn}{N})] \quad &n=0,1, \cdots, N-1 \tag{2.2} \end{aligned} \right.

x(n)x(n)XF(k)X^F(k) 均为长度为 NN 的序列。

注意k=0N1WNk=0\sum_{k=0}^{N-1} W_{N}^{k}=0
一般情况下,当 pp 为整数时,k=0N1WNpk=Nδ(p)\sum_{k=0}^{N-1} W_{N}^{p k}=N \delta(p)

# 矩阵形式

DFT:

[XF(0)XF(1)XF(k)XF(N1)]=[WN0WN0WN0WN0WN0WN1WN2WNN1WN0WN2WN4WN2(N1)WN0WNkWN2kWNk(N1)WN0WN(N1)WN2(N1)WN(N1)(N1)][x(0)x(1)x(k)x(N1)]\begin{array}{c} \left[ \begin{matrix} X^F(0) \\ X^F(1) \\ \vdots \\ X^F(k) \\ \vdots \\ X^F(N-1) \\ \end{matrix} \right] = \left[ \begin{matrix} W_{N}^{0} & W_{N}^{0} & W_{N}^{0} & \cdots & W_{N}^{0} \\ W_{N}^{0} & W_{N}^{1} & W_{N}^{2} & \cdots & W_{N}^{N-1} \\ W_{N}^{0} & W_{N}^{2} & W_{N}^{4} & \cdots & W_{N}^{2(N-1)} \\ \vdots & \vdots & \vdots & \cdots & \vdots \\ W_{N}^{0} & W_{N}^{k} & W_{N}^{2 k} & \cdots & W_{N}^{k(N-1)} \\ \vdots & \vdots & \vdots & \cdots & \vdots \\ W_{N}^{0} & W_{N}^{(N-1)} & W_{N}^{2(N-1)}& \cdots & W_{N}^{(N-1)(N-1)} \end{matrix} \right] \left[ \begin{matrix} x(0) \\ x(1) \\ \vdots \\ x(k) \\ \vdots \\ x(N-1) \\ \end{matrix} \right] \end{array}

IDFT:

[x(0)x(1)x(k)x(N1)]=1N[WN0WN0WN0WN0WN0WN1WN2WN(N1)WN0WN2WN4WN2(N1)WN0WNkWN2kWNk(N1)WN0WN(N1)WN2(N1)WN(N1)(N1)][XF(0)XF(1)XF(k)XF(N1)]\begin{array}{c} \left[ \begin{matrix} x(0) \\ x(1) \\ \vdots \\ x(k) \\ \vdots \\ x(N-1) \\ \end{matrix} \right] = \frac{1}{N} \left[ \begin{matrix} W_{N}^{0} & W_{N}^{0} & W_{N}^{0} & \cdots & W_{N}^{0} \\ W_{N}^{0} & W_{N}^{-1} & W_{N}^{-2} & \cdots & W_{N}^{-(N-1)} \\ W_{N}^{0} & W_{N}^{-2} & W_{N}^{-4} & \cdots & W_{N}^{-2(N-1)} \\ \vdots & \vdots & \vdots & \cdots & \vdots \\ W_{N}^{0} & W_{N}^{-k} & W_{N}^{-2 k} & \cdots & W_{N}^{-k(N-1)} \\ \vdots & \vdots & \vdots & \cdots & \vdots \\ W_{N}^{0} & W_{N}^{-(N-1)} & W_{N}^{-2(N-1)}& \cdots & W_{N}^{-(N-1)(N-1)} \end{matrix} \right] \left[ \begin{matrix} X^F(0) \\ X^F(1) \\ \vdots \\ X^F(k) \\ \vdots \\ X^F(N-1) \\ \end{matrix} \right] \end{array}

有以下结论:

  • DFT 矩阵 [F][F] 是对称的,即 [F]=[F]T[F]=[F]^T
  • [F][F] 是酉函数,[F][F]=N[IN][F][F]^*=N[I_N],其中 [IN][I_N](N×N)(N\times N) 的单位矩阵。

[F]1=1N[F],[F][F]1=[IN]=[1O1O1][F]^{-1}=\frac{1}{N}[F]^*,\quad [F][F]^{-1}=[I_N]= \left[ \begin{matrix} 1 & & & O \\ & 1 & & \\ & & \ddots & \\ O & & & 1 \end{matrix} \right]

于是乎,可利用矩阵化简,提高计算效率。

# 代码实现(Python)

式 (2.2) 的实现:

import numpy as np

# 一维朴素 DFT 实现
def myDFT(x):
    N = x.shape[0]
    y = np.zeros(N, dtype=complex) # 初始化
    for i in range(N):
        for k in range(N):
            real = np.cos(2 * np.pi / N * i * k)        # 实部
            imag = -1j * np.sin(2 * np.pi / N * i * k)  # 虚部
            y[i] += x[k] * (real + imag)
    return y

# 一维朴素 IDFT 实现
def myIDFT(x):
    dim = x.shape[0]
    y = np.zeros(dim, dtype=complex) # 初始化
    for i in range(dim):
        for k in range(dim):
            real = np.cos(2 * np.pi / dim * i * k)        # 实部
            imag = +1j * np.sin(2 * np.pi / dim * i * k)  # 虚部
            y[i] += x[k] * (real + imag)
        y[i] = y[i] / dim
    return y

x = np.array([1+0j,2+0j,3+0j])

print("x:      ",x)
print("np.fft: ",np.fft.fft(x))
print("myDFT:  ",myDFT(x))
print("myIDFT: ",myIDFT(myDFT(x)))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

结果:(事实上限于精度,numpy.fft.fft 与本文 myDFT 代码所得并不完全一致,可看到精度误差大概在 e16e^{-16} 级别)

>>> python DFT.py
x:       [1.+0.j 2.+0.j 3.+0.j]
np.fft:  [ 6. +0.j        -1.5+0.8660254j -1.5-0.8660254j]
myDFT:   [ 6. +0.j        -1.5+0.8660254j -1.5-0.8660254j]
myIDFT:  [1.-9.62193288e-16j 2.-1.48029737e-16j 3.+3.70074342e-16j]
1
2
3
4
5

# 推导

有关傅里叶级数推导先看看 傅里叶级数推导过程

首先对比一下连续和离散方程的差异,并注意这里都是系数

{Ak=1TTf(t)ejkωtdtXF(k)=1Nn=0N1x(n)ejkn2πN,k=0,1,2,,N1\left\{ \begin{aligned} & A_k = \frac{1}{T} \intop_{T} f(t) e^{-j k \omega t} d t \\ & X^F(k)= \frac{1}{N} \sum_{n=0}^{N-1} x(n) e^{-jkn\frac{2\pi}{N}}, \qquad k=0,1,2,\cdots,\ N-1 \end{aligned} \right.

取步长 dtTNdt \approx \frac{T}{N} 并利用梯形积分法来近似代替,并注意 ω=2πT\omega=\frac{2\pi}{T}
即有:

Ak=1T0Tf(t)ejkωtdt1TTNn=0N1f(nTN)ejk2πTnTN=1Nn=0N1x(n)ejkn2πN\begin{aligned} A_k & = \frac{1}{T} \int_{0}^{T} f(t) e^{-j k \omega t} dt \\ & \approx \frac{1}{T} \frac{T}{N} \sum_{n=0}^{N-1} f(n \cdot \frac{T}{N}) e^{-j k \frac{2\pi}{T} \ n\frac{T}{N}} \\ & = \frac{1}{N} \sum_{n=0}^{N-1} x(n) e^{-jkn \frac{2\pi}{N}} \end{aligned}

其中 x(n)=f(nTN)x(n)=f(n \cdot \frac{T}{N})

# 性质

# Z 变换

fs=1Tf_s = \dfrac{1}{T} 为采样速率,或者 T=1fsT=\dfrac{1}{f_s} 即采样间隔(单位为 s),有

X(ejωT)=n=0N1x(n)ej2πfnfsX\left(\mathrm{e}^{\mathrm{j} \omega T}\right)=\sum_{n=0}^{N-1} x(n) e^{ \frac{-\mathrm{j} 2 \pi f n}{f_{\mathrm{s}}}}

此式表示 X(z)X(z) 在 z 平面以原点为中心的单位圆上的取值。 通过对单位圆进行等间隔采样,

XF(k)=n=0N1x(n)ej2πknNk=0,1,,N1X^{\mathrm{F}}(k)=\sum_{n=0}^{N-1} x(n) e^{\frac{-\mathrm{j} 2 \pi k n}{N}} \quad k=0,1, \cdots, N-1

式中 k=Nffsk=\dfrac{Nf}{f_s}

x(n)x(n) 的 DFT 实际上是 x(n)x(n) 的 Z 变换在单位圆的等间隔采样。 XF(k)X^F(k) 表示 {x(n)}\{x(n)\} 在频率为 f=kfsNf=\dfrac{kf_s}{N} 时的 DFT,k=0,1,,N1k=0,1,\cdots,N-1。 由于存在频率混叠,因此 x(n)x(n) 的最高频率系数为 XF(N2)X^F(\dfrac{N}{2}),即 f=fs2f=\dfrac{f_s}{2}

关于频率混叠,推荐看看 这篇文章 (opens new window)

# 线性

[a1x1(n)+a2x2(n)][a1X1F(k)+a2X2F(k)]\left[a_{1} x_{1}(n)+a_{2} x_{2}(n)\right] \Leftrightarrow\left[a_{1} X_{1}^{\mathrm{F}}(k)+a_{2} X_{2}^{\mathrm{F}}(k)\right]

其中,a1,a2a_1,a_2 均为常数

# 复共轭定理

x(n)x(n) 为实序列时,有

XF(N2+k)=XF(N2k)k=0,1,,N2X^{\mathrm{F}}\left(\frac{N}{2}+k\right)=X^{\mathrm{F}*} \left(\frac{N}{2}-k\right) \quad k=0,1, \cdots, \frac{N}{2}

功率谱在频域内是关于 N2\frac{N}{2} 对称的偶函数

# 帕斯瓦尔定理

该定理对所有酉变换都成立

n=0N1x(n)x(n)=1Nk=0N1XF(k)XF(k)\sum_{n=0}^{N-1} x(n) x^{*}(n)=\frac{1}{N} \sum_{k=0}^{N-1} X^{\mathrm{F}}(k) X^{\mathrm{F}^{*}}(k)

证明:

n=0N1x(n)2=n=0N1x(n)x(n)=1Nn=0N1x(n)k=0N1XF(k)WNnk=1Nk=0N1XF(k)n=0N1x(n)WNnk=1Nk=0N1XF(k)(n=0N1x(n)WNnk)=1Nk=0N1XF(k)XF(k)=1Nk=0N1XF(k)2\begin{aligned} \sum_{n=0}^{N-1}|x(n)|^{2} &=\sum_{n=0}^{N-1} x(n) x^{*}(n)=\frac{1}{N} \sum_{n=0}^{N-1} x^{*}(n) \sum_{k=0}^{N-1} X^{\mathrm{F}}(k) W_{N}^{-n k} \\ &=\frac{1}{N} \sum_{k=0}^{N-1} X^{\mathrm{F}}(k) \sum_{n=0}^{N-1} x^{*}(n) W_{N}^{-n k}=\frac{1}{N} \sum_{k=0}^{N-1} X^{\mathrm{F}}(k)\left(\sum_{n=0}^{N-1} x(n) W_{N}^{n k}\right)^{*} \\ &=\frac{1}{N} \sum_{k=0}^{N-1} X^{\mathrm{F}}(k) X^{\mathrm{F}*}(k)=\frac{1}{N} \sum_{k=0}^{N-1}\left|X^{\mathrm{F}}(k)\right|^{2} \end{aligned}

# 循环移位

给定 x(n)XF(k)WNhkx(n)\Leftrightarrow X^{\mathrm{F}}(k)W_N^{-hk},则

x(n+h)XF(k)WNhkx(n+h) \Leftrightarrow X^{\mathrm{F}}(k) W_{N}^{-h k}

式中,x(n+h)x(n+h)便是在时域向左循环移位 hh 个采样间隔,即 {xn+h}\{x_{n+h}\}(xh,xh+1,xh+2,,xN1,x0,x1,,xh1)(x_h,x_{h+1},x_{h+2},\cdots,x_{N-1},x_0,x_{1},\cdots,x_{h-1})

因为 WNhk=1|W_N^{-hk}| = 1,所以原序列的 DFT 和循环移位后的序列的 DFT 仅相位是相关的。对一个序列进行循环移位操作,其幅度谱和功率谱均保持不变。

# 置换序列的 DFT

{x(pn)}{XF(qk)}0p,qN1\{x(p_n)\} \Longleftrightarrow \{X^F(q_k)\} \qquad 0 \leq p,q \leq N-1

其中 pnp_npnp \cdot nNN 求余的余数。同时 ppNN 互质,为整数。
其中 qkq_kqkq \cdot kNN 求余的余数。同时 qqNN 互质,为整数。

即:

AF(k)=XF(qk)=n=0N1x(pn)WnkA^{\mathrm{F}}(k)=X^F(q_k)=\sum_{n=0}^{N-1} x(p_n) W^{n k}

附一个示例:

N=8p=3N=8、p=3,则有 q=3q=3。令 x(n)=0,1,2,3,4,5,6,7x(n)={0,1,2,3,4,5,6,7},则

AF(K)=XF(qk)={XF(0),XF(3),XF(6),XF(1),XF(4),XF(7),XF(2),XF(5)}A^{\mathrm{F}}(K)=X^{\mathrm{F}}(q k)=\left\{X^{\mathrm{F}}(0), X^{\mathrm{F}}(3), X^{\mathrm{F}}(6), X^{\mathrm{F}}(1), X^{\mathrm{F}}(4), X^{\mathrm{F}}(7), X^{\mathrm{F}}(2), X^{\mathrm{F}}(5)\right\}

则有:

a(n)=[AF(k)]NIDFT=0,3,6,1,4,7,2,5=x(pn)\begin{aligned} a(n)&=[A^F(k)] 的 N 点 IDFT \\ &={0, 3, 6, 1, 4, 7, 2, 5} \\ &=x(p_n) \end{aligned}

p=5p=5 时,q=5q=5;当 p=7p=7 时,q=7q=7,同理。

本文为《快速傅里叶变换 算法与应用》的笔记,部分内容为摘抄,请悉知。

Last Updated: 2023-10-29T08:26:04.000Z