Engee 文档
Notebook

DSP基础：快速傅立叶变换（FFT）

快速傅立叶变换(FFT，FFT-Fast Fourier Transform)是对离散傅立叶变换(dft)进行快速计算的算法，广泛应用于数字信号处理(DSP)中进行频域分析处理。 在DSP中，经常需要将信号从时域转移到频域（反之亦然），以便:

*分析信号频谱
*过滤信号
*压缩数据
*实现调制算法（WI-Fi中的OFDM，5G）
*处理雷达和医疗信号

让我们来看看如何从插件库中使用FFT函数。 FFTW.jl 用于信号频谱分析的任务:

In [ ]: 
using FFTW

一个简单的信号及其频谱

让我们熟悉使用该函数的功能 fft 使用持续时间为0.2秒的真实正弦信号的示例来分析频谱，采样频率为2kHz，振幅为0.5，重复频率为100Hz。:

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, " Гц")
df = fs/nsamp;
freq_vec = -fs/2:df:fs/2-df;
plot(freq_vec, abs.(B))

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

现在我们看到频率为-100Hz和+100Hz的两个峰值。 但它们的振幅与所讨论的正弦曲线的振幅不对应。

为了准确估计信号在其频谱上的能量，我们将FFT的输出矢量归一化为样本数，即除以 nsamp. 但是，由于我们实际信号的一半能量已经"进入"负频率范围，我们需要将得到的频谱乘以2。 那么，我们将使用条形图来显示（stem)在正频率范围内的圆形标记，由0至400赫兹:

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]:

我们观察到频率为100Hz的峰值，幅度为0.5，这与原始正弦信号的参数完全对应。

更改频率分辨率

现在让我们尝试改变输入信号的采样数（FFT的"窗口"）。 让我们将函数应用于输入 fft 正弦波的前80个样本。 在这种情况下，FFT算法的输出也将具有80个复数样本，这意味着频率网格的步长将发生变化。 信号变短了5倍，因此，频率分辨率变成了25赫兹。:

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]:

改变频率网格不影响信号幅度和频率的估计。 25Hz的折扣使我们能够"准确地击中"100Hz的频谱分量。 但是，如果频谱的音高不是估计正弦波频率的倍数，会发生什么？

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]:

我们观察到频率为91Hz和109Hz的两个峰值。 这是频谱泄漏的现象：信号在100Hz时的能量被频率网格的相邻样本上的窗函数的主瓣和旁瓣（参见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.1，而不是0.5。 这是由于这样的事实，即信号已经成为五倍长，由于零，其能量，分别也被削弱了五倍。

当FFT点数为2的幂时，经典算法的工作效率更高，并且添加零对于长度对齐和可视化很有用，但它不会增加真实的频率分辨率。 对于频谱分析的真正改进，有必要增加信号的记录时间。

交互式频谱分析

最后，我们将对频率为80Hz、120Hz和245Hz，振幅分别为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 来创建代码单元掩码。 我们可以改变频率分辨率值（RBW）使用该工具 slider.

在单元代码中，计算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到20Hz，并查看估计信号的单个频谱分量的频率和幅度的准确性。

结论

我们熟悉了使用该函数的特点 fft 图书馆 FFTW.jl 来评估离散信号的频谱，并且还学习了如何应用码元的最大值并交互地执行数字信号处理任务。