快速傅立叶变换

Jia-Yin2024年3月7日大约 5 分钟

离散傅立叶变换 (Discrete Fourier Transform, DFT)

前面提到， $s(t)$ 的傅立叶变换为

S(f) = \int_{-\infty}^{\infty} s(t) e^{-j 2 \pi f t} \, dt

在许多应用当中，我们很难找到 $s(t)$ 的函数形式，即便能够找到，在计算积分的时候，也很可能找不到封闭的积分形式，因此大多数的情况之下，我们会使用数值的方式来计算傅立叶变换，最常见的计算方法，可能是离散傅立叶变换 (Discrete Fourier Transform, DFT)。

假设信号 $s(t)$ 在一个区间采样 $N$ 个点，其值为 $s_0, s_1, \cdots, s_{N-1}$ ，定义其 DFT 为：

S_k = \sum_{n=0}^{N-1} s_n e^{-j\frac{2\pi}{N}nk}

注意在上式中，当 $k$ 增加 $N$ 时，其值不变，因此主要的变换值为 $S_0, S_1, \cdots, S_{N-1}$ 。

比较上述的积分式以及 DFT 的定义，会发现两个公式有一些相似的地方，由于 DFT 处理的是采样之后的信号值，而且是只局限在一个区间的 $N$ 个点，因此我们可以把 DFT 算出来的 $S_k$ 值当作是 $S(f)$ 的近似采样值。

根据前面的采样定理，假设采样速度为 $f_s$ ，则可以还原的信号频率介于 [- $f_s/2$ , $f_s/2$ ] 之间，将此频率范围均匀分成 $N$ 个等分，每等分大小为

\Delta f = \frac{f_s/2 - (-f_s/2)}{N} = \frac{f_s}{N} = \frac1{NT}

这边 $T$ 为采样时间间隔，因此 $NT$ 为整个采样区间的总共时间。归纳以上提到的 DFT 的几个重要性质如下：

将 $N$ 个时间采样点 $s_0, s_1, \cdots, s_{N-1}$ ，变换为 $N$ 个频率采样点 $S_0, S_1, \cdots, S_{N-1}$ ，其中 $S_k$ 是周期性的，周期为 $N$ 。
假设采样速度为 $f_s$ ，则可以决定的最高频率值为 $f_s/2$ 。
采样频率之间的间隔为采样总时间长度的倒数。

DFT 还有很多其他重要的性质，这边再列出三个重要性质：

DFT 运算是线性的，亦即 DFT( $ax+by$ ) = $a$ DFT( $x$ ) + $b$ DFT( $y$ )。
如果 $s_n$ 是实数，则 $S_k = S^*_{-k} = S^*_{N-k}$ ，此处 * 表示共轭复数。
根据上述的 DFT 公式，可以得到 DFT 的反变换如下：

s_n = \sum_{k=0}^{N-1} S_k e^{j\frac{2\pi}{N}nk}

快速傅立叶变换 (Fast Fourier Transform, FFT)

在计算离散傅立叶变换的时候，可以利用公式中许多对称和共轭的性质来加快计算过程，这种快速计算的方法称为快速傅立叶变换，简称 FFT。FFT 有不同的计算形式，但实际上就是快速计算 DFT 的算法。FFT 的逆变换称为 Inverse FFT，或称为 IFFT。

一般在信号处理时，我们多使用快速傅立叶变换 (FFT) 来计算其频谱，使用快速傅立叶变换时，应牢记以下几个重点：

假设用来分析的信号段有 $N$ 个点，且其采样间隔为 $T$ ，则此信号段总长为 $NT$ ，且采样速度为 $f_s=1/T$ 。
根据采样定理，可以观察的频率范围为 $[-f_M,f_M]$ ， $f_M=f_s/2$ 。
$N$ 个时域点做 FFT 之后变成 $N$ 个频域点，故频率解析度，亦即相邻频率点的间隔为 $2f_M/N = f_s/N = 1/(NT)$ 。简单说，采样时间总长的倒数就是频率采样点之间的距离，也就是频率解析度。
要记得采样速度的一半，是可观察的最大频率；采样总长的倒数是频率解析度。
如果是实函数的话，做完 FFT 之后，其正频率和负频率相对的值会是共轭复数。

我们现在来观察帽子函数的频谱，帽子函数 rect(t) 图形如下：

使用 MATLAB/Octave 观察其频谱步骤如下：

假设采样范围从-1.5到1.5，每隔0.1取一点，先计算时域信号：

t1 = -1.5:0.1:-0.6;
t2 = -0.4:0.1:0.4;
t3 =  0.6:0.1:1.5;
t = [t1 -0.5 t2 0.5 t3];
s = [zeros(size(t1)) 0.5 ones(size(t2)) 0.5 zeros(size(t3))];
plot(t, s);

这个程式片段写得有点繁琐，不过很容易理解，执行看看结果是否正确。

计算频谱，直接使用 fft 函数将上面的时域信号转到频域，应注意转出来的结果为复数，我们先观察其振幅。

sf = fft(s);
plot(abs(sf));

注意在上面的程式中，画图时只给了y轴，横轴会是点数，得到的结果如下：

上面提过实函数的 FFT，其正频率和负频率相对的值是共轭复数。实际上我们现在看到的图形中，图的右半边应该移到负的那边去才对。MATLAB 里面有一个 fftshift 就做这事。另外，最大频率为采样率一半=5，间隔为总长倒数=1/3，我们要把横轴 (频率) 加上去。

sf = fft(s);
sf = abs(fftshift(sf));
f = -5:1/3:5;
plot(f, sf);
grid on;

最后得到修正的图如下，这个图形就是 sinc 函数的振幅。

完整程式码如下：

t1 = -1.5:0.1:-0.6;
t2 = -0.4:0.1:0.4;
t3 =  0.6:0.1:1.5;
t = [t1 -0.5 t2 0.5 t3];
s = [zeros(size(t1)) 0.5 ones(size(t2)) 0.5 zeros(size(t3))];
sf = fft(s);
sf = abs(fftshift(sf));
f = -5:1/3:5;
plot(f, sf);
grid on;

练习 4

上面的频谱图看起来并不是非常平滑，请修改程式的参数，让画出来的频谱变得更加平滑。