跳至主要內容

快速傅立叶变换

Jia-Yin大约 5 分钟

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

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

S(f)=∫−∞∞s(t)e−j2πft dt S(f) = \int_{-\infty}^{\infty} s(t) e^{-j 2 \pi f t} \, dt

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

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

Sk=∑n=0N−1sne−j2πNnk S_k = \sum_{n=0}^{N-1} s_n e^{-j\frac{2\pi}{N}nk}

注意在上式中，当 kk 增加 NN 时，其值不变，因此主要的变换值为 S0,S1,⋯ ,SN−1S_0, S_1, \cdots, S_{N-1}。

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

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

Δf=fs/2−(−fs/2)N=fsN=1NT \Delta f = \frac{f_s/2 - (-f_s/2)}{N} = \frac{f_s}{N} = \frac1{NT}

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

  1. 将 NN 个时间采样点 s0,s1,⋯ ,sN−1s_0, s_1, \cdots, s_{N-1}，变换为 NN 个频率采样点 S0,S1,⋯ ,SN−1S_0, S_1, \cdots, S_{N-1}，其中 SkS_k 是周期性的，周期为 NN。
  2. 假设采样速度为 fsf_s，则可以决定的最高频率值为 fs/2f_s/2。
  3. 采样频率之间的间隔为采样总时间长度的倒数。

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

  1. DFT 运算是线性的，亦即 DFT(ax+byax+by) = aa DFT(xx) + bb DFT(yy)。
  2. 如果 sns_n 是实数，则 Sk=S−k∗=SN−k∗S_k = S^*_{-k} = S^*_{N-k}，此处 * 表示共轭复数。
  3. 根据上述的 DFT 公式，可以得到 DFT 的反变换如下：

sn=∑k=0N−1Skej2πNnk 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) 来计算其频谱，使用快速傅立叶变换时，应牢记以下几个重点：

  1. 假设用来分析的信号段有 NN 个点，且其采样间隔为 TT，则此信号段总长为 NTNT，且采样速度为 fs=1/Tf_s=1/T。

  2. 根据采样定理，可以观察的频率范围为 [−fM,fM][-f_M,f_M]，fM=fs/2f_M=f_s/2。

  3. NN 个时域点做 FFT 之后变成 NN 个频域点，故频率解析度，亦即相邻频率点的间隔为 2fM/N=fs/N=1/(NT)2f_M/N = f_s/N = 1/(NT)。 简单说，采样时间总长的倒数就是频率采样点之间的距离，也就是频率解析度。

  4. 要记得采样速度的一半，是可观察的最大频率；采样总长的倒数是频率解析度。

  5. 如果是实函数的话，做完 FFT 之后，其正频率和负频率相对的值会是共轭复数。

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

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

  1. 假设采样范围从-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);

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

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

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

  1. 上面提过实函数的 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

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

Copyright © 2024 Jia-Yin