Engee 文档
Notebook

DSP 基础知识：快速傅立叶变换 (FFT)

快速傅立叶变换 (FFT) 是一种快速离散傅立叶变换 (DFT) 算法，广泛应用于数字信号处理 (DSP) 中的频域分析和处理。在 DSP 中，经常需要将信号从时域转换到频域（反之亦然），以便：

  • 分析信号的频谱
  • 滤波信号
  • 压缩数据
  • 执行调制算法（Wi-Fi、5G 中的 OFDM）
  • 处理雷达和医疗信号

让我们看看如何将FFTW.jl 插件库中的 FFT 函数用于信号频谱分析任务：

In [ ]: 
using FFTW

一个简单信号及其频谱

让我们以持续时间为 0.2 秒、采样率为 2 kHz、振幅为 0.5、重复频率为 100 Hz 的实际正弦信号为例，了解使用fft 功能分析频谱的特殊性：

In [ ]: 
fs = 2000;
dt = 1/fs;
t = 0:dt:0.2-dt;
sinewave = 0.5 * sin.(2*pi*t*100);
plot(t,sinewave)
Out[0]:

值得注意的是离散信号的长度sinewave - 400 个采样点。这对进一步分析非常重要。

In [ ]: 
nsamp = length(t)
Out[0]:
400

让我们应用fft 函数，将结果可视化：

In [ ]: 
A = fft(sinewave);
plot(A)
Out[0]:

FFT 函数的输出是一个复数矢量A ，其采样数目与时域中的原始信号相同。

要估算信号的振幅频谱，必须将 FFT 输出的采样取模：

In [ ]: 
plot(abs.(A))
Out[0]:

我们可以观察到对应于正负频率区域的两个峰值。但这种表示法并不能为我们提供有关这些峰值的频率或振幅的可靠信息。

为了更清晰地表示频谱，我们将使用函数fftshift 来移动 FFT 的输出矢量，同时设置一个频率矢量，并将其沿 X 轴延展。我们在所谓的第一奈奎斯特区（从-fs/2 到 +fs/2）观察信号频谱。我们可以用整个频率范围（从-fs/2 到 +fs/2，或从 0 到 fs）除以 FFT 点数，从而得出 RBW（分辨率带宽）或频率网格步长：

In [ ]: 
B = fftshift(A);
println("Разрешение по частоте (RBW) = ", fs/nsamp, " Гц")
freq_vec = -fs/2:fs/nsamp:fs/2-df;
plot(freq_vec, abs.(B))

Разрешение по частоте (RBW) = 5.0 Гц
Out[0]:

现在我们看到两个峰值，分别位于 -100 Hz 和 +100 Hz。但它们的振幅与正弦波的振幅不一致。

为了准确估计信号频谱上的能量，我们将 FFT 输出矢量按采样数归一化，即除以nsamp 。但由于实际信号的一半能量 "跑 "到了负频率区域，因此我们需要将得到的频谱乘以 2。为了便于显示，我们将使用线形图 (stem)，在 0 至 400 Hz 的正频率区域使用圆形标记：

In [ ]: 
S = 2*B/nsamp;
println("Разрешение по частоте (RBW) = ", fs/nsamp, " Гц")
plot(freq_vec, abs.(S), 
        line=:stem,
        marker=:circle,
        xlim = (0, 400),
        title = "Амплитудный спектр синусоиды")

Разрешение по частоте (RBW) = 5.0 Гц
Out[0]:

我们在 100 Hz 处观察到一个振幅为 0.5 的峰值，与原始正弦信号的参数完全一致。

频率分辨率的变化

现在，让我们尝试改变输入信号的样本数（FFT 的 "窗口"）。让我们向函数fft 输入正弦波的前 80 个采样点。在这种情况下，FFT 算法的输出也将是 80 个复数样本，这意味着频率网格的间距将发生变化。信号变短了 5 倍，因此频率分辨率变为 25 Hz：

In [ ]: 
winsize = 80;
rbw80 = fs/winsize;
freq80 = 0:rbw80:fs - rbw80;
S80 = (2*fft(sinewave[1:winsize]))./winsize;
println("Разрешение по частоте (RBW) = ", rbw80, " Гц")
plot(freq80, abs.(S80), 
        line=:stem,
        marker=:circle,
        xlim = (0, 400))

Разрешение по частоте (RBW) = 25.0 Гц
Out[0]:

频率网格的变化不会影响对信号振幅和频率的估计。25 赫兹的离散度使我们能够 "精确地 "捕捉到 100 赫兹的频谱成分。但是，如果频谱阶跃不是估计正弦波频率的倍数，会发生什么情况呢？

In [ ]: 
winsize = 110;
rbw110 = fs/winsize;
freq110 = 0:rbw110:fs - rbw110;
S110 = (2*fft(sinewave[1:winsize]))./winsize;
println("Разрешение по частоте (RBW) = ", rbw110, " Гц")
plot(freq110, abs.(S110), 
        line=:stem,
        marker=:circle,
        xlim = (0, 400))

Разрешение по частоте (RBW) = 18.181818181818183 Гц
Out[0]:

我们在 91 赫兹和 109 赫兹处观察到两个峰值。这是一种频谱泄漏现象：100 Hz 处的信号能量已被频率网格相邻样本上窗口函数的主叶和边叶（见 DSP 理论）捕获。

在 FFT 中添加零点

FFT 中的 "零填充 "技术是指在执行 FFT 之前，将原始信号填充为所需长度的零。添加零点可以增加 FFT 点的数量，降低频率网格的间距，使频谱在视觉上更加平滑。但这并不会增加信号的新信息，而只是对现有频谱进行插值。

让我们在正弦波中加入零，直到矢量长度达到 2000 个采样点：

In [ ]: 
winsize = 2000;
rbw2000 = fs/winsize;
freq2000 = 0:rbw2000:fs - rbw2000;
numofzeros = winsize - length(sinewave);
println("Добавлено нулей = ", numofzeros)
S2000 = (2*fft([sinewave; zeros(numofzeros)]))./winsize;
println("Разрешение по частоте (RBW) = ", rbw2000, " Гц")
plot(freq2000, abs.(S2000), 
        line=:stem,
        marker=:circle,
        xlim = (0, 400))

Добавлено нулей = 1600
Разрешение по частоте (RBW) = 1.0 Гц
Out[0]:

我们可以发现，频谱的峰值从 0.5 变成了 0.1。这是由于零点使信号长度增加了五倍，其能量也相应衰减了五倍。

当 FFT 点数为 2 时，经典算法的工作效率更高，增加零点有助于长度均衡和可视化，但并不能提高真正的频率分辨率。要真正提高频谱分析能力，需要延长信号记录的时间。

交互式频谱分析

最后，让我们将频谱估算技术应用于一个信号，它是频率分别为 80 Hz、120 Hz 和 245 Hz，振幅分别为 10、5 和 2.4 的三个正弦波之和。让我们在时域中显示所得矢量的形状：

In [ ]: 
three_sines = [10 5 2.4] .* sin.(2*pi*t*[80 120 245]);
sig = sum(three_sines, dims = 2);
sig = vec(sig);
plot(t,sig)
Out[0]:

现在让我们使用Engee 的脚本功能来创建代码单元掩码。我们可以使用工具slider 更改频率分辨率的值 (RBW) 。

代码单元计算 FFT 输出点数，将其与离散信号在时间上的点数进行比较，如果一个点数超过另 一个点数，则执行加零或截断算法的输入 "窗口"。然后计算并绘制振幅频谱。

In [ ]: 
RBW = 10 # @param {type:"slider",min:1,max:20,step:1}
nfft = Int64(round(fs/RBW));
println("Точек БПФ:  ", nfft,
"    Отсчётов сигнала:  ", nsamp);

if nsamp >= nfft
        win = sig[1:nfft];
        Y = 2*fft(win)./nfft;
        println("Усечение входного окна")
else
        win = [sig; zeros(nfft - nsamp)];
        Y = 2*fft(win)./nsamp;
        println("Дополнение нулями: интерполяция спектра")
end

ff = (0:nfft-1)*fs/nfft;
plot(ff, abs.(Y), 
        line=:stem, 
        xlim = (0, 400), 
        marker =:circle,
        lab = "Сумма синусоид",
        title = "Амплитудный спектр сигнала")

Точек БПФ:  200    Отсчётов сигнала:  400
Усечение входного окна
Out[0]:

尝试使用RBW 的不同值（从 1 Hz 到 20 Hz），看看对信号各个频谱成分的频率和振幅估计的准确性。

结论

我们了解了使用FFTW.jl 库中的函数fft 估算离散信号频谱的特殊性，还学习了如何在交互模式下应用代码单元的最大值和执行数字信号处理任务。