傅里叶变换（DFT）从原理到代码

离散傅里叶变换（DFT）

一、傅里叶变换公式

1. 连续傅里叶变换
   连续傅里叶变换是使用一对正交基$cos\omega t$和$sin\omega t$对时域序列进行“探测”（求时域序列与正交基的相关性），筛选出与原时域序列相关的正交基，即为原时域序列的频域序列。连续傅里叶变换公式如下：
   $$F(\omega )={\int }_{-\infty }^{+\infty }f(t)e^{-i\omega t}dt$$
2. 离散傅里叶变换
   离散傅里叶变换（DFT）方法，是以频率$f_s$对输入信号对输入信号进行等间隔采样n次后，对这n个离散点进行DFT处理得到n个复数，每个复数表示以$f_s/n$为间隔的频率成分的幅度和相位信息，即获得输入信号的频谱信息。
   对于一个长度为N的时域离散序列$x(n)$，DFT变换结果为$X(k)$，则：
   X ( k ) = D F T { x ( n ) } = 1 N ∑ n = 0 N − 1 x ( n ) × e − j 2 π N n k , k = 0 , 1 , 2 , . . . , N − 1 X(k)=DFT\{ x(n)\} = \frac{1}{N}\sum\limits_{n = 0}^{N - 1} {x(n) \times {e^{ - j\frac{{2\pi }}{N}nk}}} ,k = 0,1,2,...,N - 1 X(k)=DFT{x(n)}=N1​n=0∑N−1​x(n)×e−jN2π​nk,k=0,1,2,...,N−1
   转换为矩阵形式：
   D F T { x ( n ) } = 1 N [ 1 1 1 ⋯ 1 1 e − j 2 π N × 1 e − j 2 π N × 2 ⋯ e − j 2 π N × ( N − 1 ) 1 e − j 2 π N × 2 × 1 e − j 2 π N × 2 × 2 ⋯ e − j 2 π N × 2 × ( N − 1 ) ⋮ ⋮ ⋮ ⋱ ⋮ 1 e − j 2 π N × ( N − 1 ) × 1 e − j 2 π N × ( N − 1 ) × 2 ⋯ e − j 2 π N × ( N − 1 ) 2 ] [ x ( 0 ) x ( 1 ) x ( 2 ) ⋮ x ( N − 1 ) ] DFT\{ x(n) \}=\frac{1}{N} \begin{bmatrix} 1&1&1& \cdots &1\\ 1&{{{\rm{e}}^{ - j\frac{{2\pi }}{N} \times 1}}}&{{{\rm{e}}^{ - j\frac{{2\pi }}{N} \times 2}}}& \cdots &{{{\rm{e}}^{ - j\frac{{2\pi }}{N} \times (N - 1)}}}\\ 1&{{{\rm{e}}^{ - j\frac{{2\pi }}{N} \times 2 \times 1}}}&{{{\rm{e}}^{ - j\frac{{2\pi }}{N} \times 2 \times 2}}}& \cdots &{{{\rm{e}}^{ - j\frac{{2\pi }}{N} \times 2 \times (N - 1)}}}\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1&{{{\rm{e}}^{ - j\frac{{2\pi }}{N} \times (N - 1) \times 1}}}&{{{\rm{e}}^{ - j\frac{{2\pi }}{N} \times (N - 1) \times 2}}}& \cdots &{{{\rm{e}}^{ - j\frac{{2\pi }}{N} \times {{(N - 1)}^2}}}} \end{bmatrix} \begin{bmatrix} {x(0)}\\ {x(1)}\\ {x(2)}\\ \vdots \\ {x(N - 1)} \end{bmatrix} DFT{x(n)}=N1​⎣⎢⎢⎢⎢⎢⎡​111⋮1​1e−jN2π​×1e−jN2π​×2×1⋮e−jN2π​×(N−1)×1​1e−jN2π​×2e−jN2π​×2×2⋮e−jN2π​×(N−1)×2​⋯⋯⋯⋱⋯​1e−jN2π​×(N−1)e−jN2π​×2×(N−1)⋮e−jN2π​×(N−1)2​⎦⎥⎥⎥⎥⎥⎤​⎣⎢⎢⎢⎢⎢⎡​x(0)x(1)x(2)⋮x(N−1)​⎦⎥⎥⎥⎥⎥⎤​

二、离散傅里叶变换实践

设原信号为$f(t)=\cos (2\pi t+0.5)+3.2\sin (\pi t+1)+0.3$，时间跨度为0到$T=10s$，信号中最高频率为1Hz，则使用3Hz（两倍最高频率）对原信号进行采样，则：
注：角频率 $\omega =2\pi f$，$f$为频率。

import matplotlib.pyplot as plt
import numpy as np

fs = 30  # 采样频率为3Hz，可变更，但需要大于2（采样定理）
T = 10  # 时序样本长度为10s
N = T * fs  # 采样次数
t = np.arange(0, T, 1 / fs)
tim = np.cos(2 * np.pi * t + 0.5) + 3.2 * np.sin(np.pi * t + 1) + 0.3  # 生成时序样本

X = np.dot(np.arange(N).reshape((-1, 1)), np.arange(N).reshape((1, -1)))
X = np.exp(-1j * 2 * np.pi / N * X)
freq = np.dot(X, tim)

omega = np.arange(0, fs, fs / N)  # 生成频域图横坐标，单位Hz

filtering = np.where(abs(freq) > 0.000001, 1, 0)  # 过滤器，过滤数值奇异的零点
amplitude = abs(freq) / N * filtering  # 此处的幅度，直流项和非直流项存在差异，直流项为原值，非直流项为原幅值的一半
phase = np.angle(freq) * filtering

fig1, axes = plt.subplots(2, 1)
axes[0].bar(omega, height=amplitude, width=fs / N, label='amplitude')  # 幅值图
axes[0].legend()

axes[1].bar(omega, height=phase, width=fs / N, label='phase')  # 相位图
axes[1].legend()
plt.show()

fig2, axes = plt.subplots(4)
axes[0].plot(t, np.cos(2 * np.pi * t + 0.5))  # 原信号第一个分量
axes[1].plot(t, 3.2 * np.sin(np.pi * t + 1))  # 原信号第二个分量
axes[2].plot(t, np.ones_like(t) * 0.3)  # 原信号直流项
axes[3].plot(t, tim)

tim_back = np.zeros_like(t)
for i in range(len(omega)):
    tim_back = tim_back + amplitude[i] * np.cos(2 * np.pi * omega[i] * t + phase[i])

axes[3].plot(t, tim_back, '--', marker='*')  # 复原原图形
plt.show()

以上为根据公式，自行推导的傅里叶变换及逆变换代码，以下采用Python自带的fft工具对以上变换进行验证：

import matplotlib.pyplot as plt
import numpy as np

fs = 30  # 采样频率为3Hz，可变更，但需要大于2（采样定理）
T = 10  # 时序样本长度为10s
N = T * fs  # 采样次数
t = np.arange(0, T, 1 / fs)
tim = np.cos(2 * np.pi * t + 0.5) + 3.2 * np.sin(np.pi * t + 1) + 0.3  # 生成时序样本

"""
以下2行为Python自带FFT程序
"""
omega = np.fft.fftfreq(N, d=1 / fs)
freq = np.fft.fft(tim)

filtering = np.where(abs(freq) > 0.000001, 1, 0)  # 过滤器，过滤数值奇异的零点
amplitude = abs(freq) / N * filtering  # 此处的幅度，直流项和非直流项存在差异，直流项为原值，非直流项为原幅值的一半
phase = np.angle(freq) * filtering

fig1, axes = plt.subplots(2, 1)
axes[0].bar(omega, height=amplitude, width=fs / N, label='amplitude')  # 幅值图
axes[0].legend()

axes[1].bar(omega, height=phase, width=fs / N, label='phase')  # 相位图
axes[1].legend()
plt.show()

fig2, axes = plt.subplots(4)
axes[0].plot(t, np.cos(2 * np.pi * t + 0.5))  # 原信号第一个分量
axes[1].plot(t, 3.2 * np.sin(np.pi * t + 1))  # 原信号第二个分量
axes[2].plot(t, np.ones_like(t) * 0.3)  # 原信号直流项
axes[3].plot(t, tim)

tim_back = np.zeros_like(t)
for i in range(len(omega)):
    tim_back = tim_back + amplitude[i] * np.cos(2 * np.pi * omega[i] * t + phase[i])

"""
以下1行为Python自带FFT程序
"""
tim_back2 = np.fft.ifft(freq)

axes[3].plot(t, tim_back, '--', marker='*')  # 复原原图形
axes[3].plot(t, tim_back2.real, '--', marker='+')  # 复原原图形
plt.show()

以上代码，相较于自制的FFT程序，不同在于，Python自带np.fft.fft工具的频率从负到正，其余一致，如下：

$$图：python自带numpy.fft工具$$

$$图：自编程序的fft工具$$
原因在于DFT公式，具有$N$的周期性，如下：
$$e^{-j\frac{2\pi }{N}n(k-N)}=e^{-j\frac{2\pi }{N}nk+j2\pi n}=e^{-j\frac{2\pi }{N}nk}{e}^{j2\pi n}=e^{-j\frac{2\pi }{N}nk}$$
又因为DFT函数具有共轭对称性，因此可以将观测周期从$0\sim N$变更到$-N/2\sim N/2$之间，结果不受影响。

三、numpy.fft应用

为简洁起见，编程中常使用numpy.fft，因此以下说明numpy.fft的使用方法。

1. 离散序列输入
   设原始时序为：
   $$f(t),t\in [0,T]$$
   采样频率为$f_s$，则采样数目$N=T\ast f_s$。
2. 傅里叶变换
   采样numpy.fft进行傅里叶变换：
   获取频率分布：$freq=numpy.fft.fftfreq(N,d=1/f_s)$
   获取傅里叶变换结果，结果为复数： r e = n u m p y . f f t . f f t { f ( t ) } re=numpy.fft.fft\{f(t) \} re=numpy.fft.fft{f(t)}
   获取幅值分布：$amplitude=abs(re)/N$
   获取初相分布：$phase=np.angle(re)$

---

注：

Ⅰ. 可以使用 $filtering=np.where(abs(freq)>0.000001,1,0)$过滤掉数值奇异的数据：
$phase=phase\ast filtering$
$amplitude=amplitude\ast filtering$
Ⅱ. 幅值分布，需要÷N，即 $amplitude=abs(re)/N$

---

3. numpy.fft傅里叶结果
   numpy.fft结果，频率分布为$-N/2\sim N/2$，序列为：
   $$freq=[0,1,...,\frac{N-1}{2},-\frac{N}{2},-\frac{N-1}{2},...,-1]$$
   Python实践中的频率分布如下：
   其中频率=0表示直流分量，其幅度为直流分量的大小；频率≠0时，幅度为相应频率下幅度的1/2。
   如下图：
   
4. numpy.fft.rfft结果
   如果只需要保留正频率部分，可以使用numpy.fft.rfft。
   rfft和fft的结果对比如下：